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Abstract. Wc consider the time evolution of order parameter correlation functions after a 
sudden quantum quench of the magnetic field in the transverse field Ising chain. Using two novel 
methods based on determinants and form factor sums respectively, we derive analytic expressions 
for the asymptotic behaviour of one and two point correlators. We discuss quenches within the 
ordered and disordered phases as well as quenches between the phases and to the quantum critical 
point. We give detailed account of both methods. 
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1. Introduction 

The non-equilibrium dynamics of isolated quantum systems after a sudden "quench" of a parameter 
characterizing their respective Hamiltonians is a subject currently under intensive theoretical and 
experimental investigation. Recent experiments on trapped ultra-cold atomic gases [1, 2, 3, 4, 5, 6] 
have established that these systems are sufficiently weakly coupled to their environments as to 
allow the observation of essentially unitary non-equilibrium time evolution on very long time 
scales. This in turn provides the opportunity of investigating fundamental questions of many-body 
quantum mechanics, which in standard condensed matter systems arc obscured by decoherence and 
dissipation. Two of the main questions raised in these works are (i) how fast correlations spread 
in quantum many-body systems and (ii) whether obscrvables such as (multi-point) correlation 
functions generically relax to time independent values, and if they do, what principles determine 
their stationary properties. The first issue was addressed in a seminal work by Lieb and Robinson [7], 
which established that in lattice many-body systems information has a finite speed of propagation 
and provided a bound on the maximal group velocity. In recent years an effective light-cone scenario 
has been proposed theoretically [8, 9, 10, 11, 12], was tested in many numerical computations 
[13, 14, 15, 16, 17, 18] and was finally observed in cold-atom experiments [4]. 

The relaxation towards stationary behaviour at first appears very surprising, because unitary 
time evolution maintains the system in a pure state at all times. The resolution of this apparent 
paradox is that in the thermodynamic limit, (finite) subsystems can and do display correlations 
characterized by a mixed state, namely the one obtained by tracing out the degrees of freedom 
outside the subsystem itself. In other words, the system acts as its own bath. 

In groundbreaking ("quantum Newton's cradle") experiments on the relaxation towards 
stationary states in systems of ultra-cold atoms Kinoshita, Wenger and Weiss [6] demonstrated 
the importance of dimensionality and conservation laws for many-body quantum dynamics out 
of equilibrium. In essence, these experiments show that three dimensional condensates relax 
quickly to a stationary state that is characterized by an effective temperature, a process known as 
"thermalization" , whereas the relaxation of quasi one-dimensional systems is slow and towards an 
unusual non-thermal distribution. This difference has been attributed to the existence of additional 
(approximate) local conservation laws in the quasi-lD case, which are argued to constrain the 
dynamics in analogy with classical integrable systems. The findings of Rcf. [6] sparked a tremendous 
theoretical effort aimed at clarifying the effects of quantum integrability on the non-equilibrium 
evolution in many-particle quantum systems [19, 20, 21, 9, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 
33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48]. Many of these studies are compatible with 
the widely held belief (see e.g. [19] for a comprehensive summary) that the reduced density matrix 
of any finite subsystem (which determines correlation functions of any local observables within the 
subsystem) of an infinite system can be described in terms of either an effective thermal (Gibbs) 
distribution or a so-called generalized Gibbs ensemble [20]. It has been conjectured that the latter 
arises for integrable models, while the former is obtained for generic systems and evidence supporting 
this view has been obtained in a number of examples [20, 21, 22, 23, 24, 25, 26, 29, 31, 32, 33]. On 
the other hand, several numerical studies [34, 36, 37, 44, 45] suggest that the full picture may well 
be more complex. 

Moreover, open questions remain even with regards to the very existence of stationary states. 
For example, the order parameter of certain mean-field models have been shown to display persistent 
oscillations [49, 50, 51, 52, 53, 54]. Non-decaying oscillations have also been observed numerically 
[36] in some non-integrable one-dimensional systems. This has given rise to the concept of "weak 
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thcrmalization" , which refers to a situation where time-averaged quantities arc thermal, while 
oscillations persist on short time-scales. We note that questions related to thcrmalization and 
spreading of correlations have also been studied using holographic techniques [55, 56, 57]. 

Given that stationary behaviour is strictly speaking a property at infinite times in the 
thermodynamic limit (and these limits do not commute) it is important to have available analytic 
results that become exact in certain limits. To that end we consider here the non-equilibrium 
dynamics of the transverse field Ising chain 

L 

H(h) = -J^2[^+i + h^], (1) 

3 = 1 

where cr" are the Pauli matrices at site j, J > and we impose periodic boundary conditions 
a L+i — a \- The Hamiltonian (1) exhibits a Z 2 symmetry of rotations around the z-axis in spin 
space by 180 degrees 

a j -> ~ a j > a = x ' V > a% 3 -> a 3- ( 2 ) 
The model (1) is a crucial paradigm of quantum critical behaviour and quantum phase 
transitions [58]. At zero temperature and in the thermodynamic limit it exhibits ferromagnetic 
(h < 1) and paramagnetic (h > 1) phases, separated by a quantum critical point at h c = 1. 
For h < 1 and L — > oo there are two degenerate ground states related by the Z 2 symmetry (2). 
Spontaneous symmetry breaking selects a unique ground state, in which spins align along the in- 
direction. On the other hand, for magnetic fields h > 1 the ground state is non-degenerate and as the 
magnetic field h is increased spins align more and more along the z-direction. The order parameter 
for the quantum phase transition is the ground state expectation value (ctJ). In the following we 
will refer to the ferromagnetic phase as the ordered phase, and to the paramagnetic one as the 
disordered phase. We note that the model (1) is (approximately) realized both in solids [59] and 
in systems of cold Rubidium atoms confined in an optical lattice [60] . In the latter realization it is 
possible to investigate its non-equilibrium dynamics experimentally. 

As shown in Appendix A, the model can be diagonalized by a Jordan- Wigner transformation 
which maps the model to spinless fermions with local annihilation operators Cj, followed by a 
Fourier transform and a Bogolioubov transformation. In terms of the momentum space Bogoliubov 
fermions a k the Hamiltonian reads 

H{h) = ^e h {k)a\a k , e h {k) = 2J^1 + h? - 2/icos(fe). (3) 

k 

For details and more precise definitions see Appendix A. 
1.1. Quench protocol and observables 

In the following we focus on a global quantum quench of the magnetic field. We assume that the 
many-body system is prepared in the ground state \^o) of Hamiltonian H{ho). At time t = 
the field hg is changed instantaneously to a different value h and one then considers the unitary 
time evolution of the system characterized by the new Hamiltonian H(h), i.e. the initial state |^o) 
evolves as 

|*o(t)>=e-" H W|* >. (4) 
The above protocol corresponds to an experimental situation, in which a system parameter has 
been changed on a time scale that is small compared to any other time scale in the system. We 
note that this can be achieved in cold-atom experiments [1, 3, 4]. 
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In this paper we study the one and two-point functions of the order parameter 

9 {t) - <* (*)|*o(t» ' ( j 

pXX ^ t] = ^mmm)^ ' pTm = ^ ~ {pX[t)f - (6) 

Although quantum quenches in the ID Ising model have been the subject of many works 
[61, 62, 63, 64, 65, 66, 67, 68, 69, 70], results on the time evolution of order parameter correlation 
functions were reported only very recently by us in a short communication [71]. In the following 
we give detailed derivations on the full asymptotic time and distance dependence of one- and two- 
point order parameter correlation functions in the thermodynamic limit for quenches within the 
two phases. These have been obtained by two novel, complementary methods, and we discuss 
both of them in detail. The first method is based on the determinant representation of correlation 
functions characteristic of free-fermionic theories. The second is based on the form-factor approach 
[72, 73, 74, 75, 76] and is applicable more generally to integrable quenches in interacting quantum 
field theories [77, 32]. We stress that this approach is qualitatively different from numerical 
approaches based on quantum integrability [78, 79, 80, 44, 46]. In particular it allows us to obtain 
analytical results in the thermodynamic limit. New results not previously reported in [71] include 
expressions for the two-point function of the order parameter for quenches within the disordered 
phase and for quenches across the critical point. 



1.2. The quench variables 

As shown in Appendix A both the initial and final Hamiltonians can be diagonalized by combined 
Jordan- Wigncr and Bogoliubov transformations with Bogoliubov angles 9® and 9^ respectively. The 
corresponding Bogoliubov fermions are related by a linear transformation (B.l) characterized by 
the difference A& = 9k — 9®. In order to parametrize the quench it is useful to introduce the quantity 

hho - (h + h ) cosfc + 1 
cosA fc = — — . (7 

y/1 + h 2 - 2/icos(fc)^l + h%- 2h cos(k) 

We note that cos is invariant under the two transformations 

(h ,h) -> (h,h ) and (ho, h) -> ^ . (8) 

However, we stress that the quantum quench itself is not invariant under the maps (8) . In the form 
factor approach a more natural quantity to consider is 

K(p) =tan[A p /2]. (9) 



1.3. Time scales for two-point correlators. 

In this manuscript we are mainly concerned with equal time correlation functions of spins separated 
by a distance £, which we take to be much larger than the lattice spacing, i.e. £ ^> 1. For fixed £ 
the time evolution is naturally divided into three regimes, which are determined by the propagation 
velocity of elementary excitations of the post-quench Hamiltonian v(k) = de ^ fc - ) . For a given final 
magnetic field h, the maximal propagation velocity is 

Wmax = max \e' h (k)\ = 2 J mm[h, 1] . (10) 
The three different regimes are: 
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Figure 1. Left panel: for intermediate times Vm^t ~ £ the behaviour of p aa (£,t) is most 
conveniently determined by considering its asyptotic expansion around infinity ("space-time 
scaling limit") along the ray D max £ = kI- This viewpoint is appropriate for any large, finite 
t and t. Right panel: the asymptotic late-time regime is reached by considering time evolution at 
fixed £. To describe this regime one should consider an asymptotic expansion of p aa (£, i) around 
t = co at fixed £. 



• Short-times v max t -C I. 

• Intermediate times u max ^ ~ £■ This regime is of particular importance for both experiments and 
numerical computations. A convenient way of describing this regime is to consider evolution 
along a particular "ray" k£ = Vm^t in space-time, see Fig. 1. In order to obtain a very accurate 
description of the dynamics at a particular point along this ray, one may then construct an 
asymptotic expansion in the single variable £ around the space-time scaling limit w max t, £ — > oo, 
k fixed. 

• Late times u ma x^ 3> This includes the limit t oo at fixed but large £. In this regime it is no 
longer convenient to consider evolution along a particular ray in space-time. In order to obtain 
accurate results for the late time dynamics, one should construct an asymptotic expansion in 
t around infinity, see Fig. 1. 

Because of the horizon effect [9, 22] the short-time regime does not display interesting features. 
We therefore focus on intermediate and late times. The most convenient way of analyzing the 
intermediate time regime is via the space-time scaling limit. It is important to note that taking 
the limit k — > of results obtained in the space-time scaling limit is not expected to reproduce 
the short-time regime. Similarly, in general taking n — > oo in the space-time scaling limit does not 
necessarily reproduce the late time behaviour. 

1.4- Organization of the manuscript. 

This is the first in a series of two papers of the dynamics in the transverse field Ising chain after a 
sudden quench of the magnetic field. The second paper, henceforth referred to as "paper II" , gives 
a detailed account of properties in the stationary state, i.e. at t = oo. The present manuscript 
deals with the time evolution of observables and is organized in the following way. In section 2 we 
present a detailed summary of our main results. In section 3 we discuss the determinant approach 
for calculating two-point functions after quantum quenches in models with free fermionic spectra. 
Section 4 introduces the form factor approach to correlation functions in integrable models after 
quantum quenches and gives a detailed account of its application to one and two point functions in 
the Ising chain. The scaling limit for quenches close to the critical point is constructed in section 
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5, and section 6 contains our conclusions. Our conventions for diagonalizing the Ising Hamiltonian 
with periodic boundary conditions are summarized in Appendix A and the initial state is expressed 
in terms of eigenstates of the post-quench Hamiltonian in Appendix B. Finally, Appendix C and 
Appendix D deal with certain technical issues. 



2. Summary and discussions of the results 

In this section we present a comprehensive summary of our main results. Our analytic results are 
valid in the thermodynamic limit L — > oo and are obtained by two different methods. 

(i) The first is based on representing correlation functions as determinants and then determining 
the asymptotic behaviour in the space-time scaling limit 

M^oo fixed. (11) 

Results obtained by this methods are exact. 

(ii) The second method employs a Lehmann representation for correlation functions, which provides 
an expansion of the correlator in powers of the functions K(k) (9). The Lehmann representation 
is recast as a low-density expansion and then the dominant terms at late times and large 
distances are summed to all orders in K(k). The method exploits the existence of a small 
parameter, namely the average densities n(k) of elementary excitations of the post-quench 
Hamiltonian H(h) with momentum k in the initial state |\&o(0)) 

(M0)Hk)\M0)) K 2 (k) 

(M>o (0)1*0(0)) 1 + K*(k) ■ [ > 

A small quench is defined as being such that these densities are small for all fc, i.e. 
maxfc_ftT 2 (/c) <C 1 (this does not necessarily imply that h and ho have to be very close to each 
other). The form factor approach provides accurate results for small quenches at late times 
and large distances. In practice the form factor method provides very good approximations to 
the exact result, except for quenches to or from the close vicinity of the critical point. 

Analytic results obtained by the two methods are compared to a direct numerical evaluation of the 
determinant representations of correlation functions in the thermodynamic limit. Finite-size effects 
are concomitantly absent and for the purposes of the comparisons shown in the following figures 
the numerical results can therefore be considered to be exact. 



2.1. One-point correlation function 



For quenches starting in the disordered phase, i.e. h > 1, the order parameter expectation value 
is zero for all times because the Z2 symmetry remains unbroken. For quenches that start and 
end in the ordered phase, i.e. h ,h < 1, exact calculations based on the determinant approach 
show that the order parameter relaxes to zero exponentially at late times > 1) (note that since 
In |cos Afe| < the exponential is always decreasing) 

dfc 

-4(fc)ln|cosA fc | , (13) 



«W)^(^ F ) 5 exp 



/ 

Jo 



7T 



where 



•-FF 



1 - fafep + V(l - - fag) 

2^1-^0(1-^)* 



(14) 
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Figure 2. Expectation value of the order parameter after a quench within the ferromagnetic 
phase. Numerical data are compared with the asymptotic predictions from both determinants and 
form factors at order 0(K 2 ). The right panel shows the accuracy of the determinant result even 
very close to the critical point, while the form factor result ceases to be an accurate approximation 
because the density of excitations is no longer small. Numerical results are obtained by considering 
the cluster decomposition of the two-point function at distance I = 180. The left hand panel 
shows that the short-time behaviour after the quench is sensitive to whether the transverse field 
is increased or decreased. If 1 > h > ho then (a x (t)) decreases for short times, while it initially 
increases if 1 > ho > h. 



This result is obtained by applying the cluster decomposition principle to the two-point function 
(19) (see Eq. (69) in Section 3 for the proof of Eq. (19)). The form factor approach gives (see 
Section 4.2) 

■)sexpU r—e' h (k) 2K 2 (k)] . (15) 



We note that 



«(()>^(1 



In | cos AJ = In 



ff) 



(1 



l-K 2 {k) 
K 2 (k) 
(h- 



1 - 

2\i 



ho) 



-2K 2 (k) - -K 6 (k) + 
3 



+ 



(16) 



64(1 -h 2 y 

so that (15) is indeed the "low-density" approximation to (13). 

Fig. 2 shows the comparison of the asymptotic result against exact numerical computation 
(obtained from cluster decomposition of the two-point function). It is evident that the asymptotic 
results become accurate already for small values of t. Eqn (13) is asymptotically valid also for 
quenches to the critical point as is shown in the left panel of Fig. 2. 

For a quench from the ferromagnetic ho < 1 to the paramagnetic h > 1 phase, we conjecture 
that the expectation value of the order parameter at late times t is given by 

i i r dk 

(af(t)) = (Cf P ) 3 [1 + cos(2e h (k )t + «) + ...] 5 exp t / — e' h {k) In |cos A fc | , (17) 



where fco is a solution of the equation cos Afc = 0, a(h, ho) is an unknown constant, and 

'hy/X 



c 



FP 



h 2 



h + h Q 



(18) 
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Figure 3. Expectation value of the order parameter after a quench from the ferromagnetic to 
the paramagnetic phase (inset) and its absolute value (main plot) in logarithmic scale to show the 
exponential decay of the correlation. Numerical data from ho = 1/2 to h = 5/3 and h = 1.05 are 
compared with the asymptotic conjecture in Eq. (17) which is reported as a continuous line. The 
asymptotic formula is valid all the way up to the critical point. Data are obtained by considering 
the cluster decomposition of the two-point function at distance i = 180. 



The dots in eqn (17) indicate subleading contributions. The conjecture (17) is compared to the 
numerically calculated one-point function (the 1-point function is obtained by applying the cluster 
decomposition principle to the two-point function Eq. (32), see Section 3.2.1) in Fig. 3. The 
agreement is clearly excellent. From a mathematical point of view the oscillating factor is a 
correction to the asymptotic behaviour, as is most clearly seen by considering log \{o~f (t)\). However, 
by virtue of its oscillatory nature, its presence obscures the leading behaviour and needs to be 
included in order to have a good description of the quench dynamics. In the limit h — > 1, goes 
to and £i(ko) — 0, signaling that the crossover between (13) and (17) is smooth. In particular, 
when approaching the critical point, the oscillation frequency decreases as shown in Fig. 3. 



2.2. Equal time two-point correlation function 



2.2.1. Quench within the ferromagnetic phase. For a quench within the ordered phase, the 

determinant approach (detailed in Section 3) leads to the following result for the two-point function 

in the space-time scaling limit (t, I — > oo with v ma ^t/H fixed) 

r 7r A}. -i 
~ \n\cos& k \0 H (2s' h {k)t-l) 



p x p F (l,t) ~ Cp F exp 



2t 



TV 

71 dfc 



n 



x exp 

Here Oh(x) is the Heaviside step function 

' 1 if x > 0, 



7T 



e' h {k)\n\cosk k \6 H (t-2e' h {k)t) 







else 



(19) 



(20) 



The constant Cp F (14) is fixed by matching (19) to the corresponding result at infinite time t — oo, 
which is derived in paper II [81]. 
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Figure 4. Numerical data for a quench within the ferromagnetic phase from ho = 1/3 to 
h = 2/3. Left: The two-point function against the asymptotic prediction Eq. (19) for i = 30 (up 
to a multiplicative factor) showing excellent agreement in the scaling regime. Inset: Ratio between 
the numerical data and asymptotic prediction (69). The leading correction is time independent, 
but subleading contributions oscillate. Right: The connected correlation function for the same 
parameters as on the left. For t < tp, p^ x (£,t) vanishes identically in the scaling regime. 



In the limit £ — > oo (19) gives the square of the result (13) for the one-point function. For 
times smaller than the Fermi time 

t F = (21) 

the first exponential factor in (19) equals 1. Thus, in the space-time scaling limit, connected 
correlations vanish identically for times t < tp and begin to form only after the Fermi time. This 
is a general feature of quantum quenches [9, 22] and has been recently observed in experiments on 
one dimensional cold-atomic gases [4]. We stress that this by no means implies that the connected 
correlations are exactly zero for t < tp: in any model, both on the lattice or in the continuum 
there are exponentially suppressed terms (in £) which vanish in the scaling limit. The form factor 
approach gives the following result for large t and I (see Section 4.3) 

pp x F (e, t) ~ (1 - h 2 ) 1 * cxp \-2t f —K 2 (k)9 H (2e , h (k)t - i) 

xexp -it —e , h (k)K 2 (k)9 H (£-2s' h (k)t). (22) 

As expected, it gives the low density approximation to the full result (19). 

A comparison (for a typical quench from ho = 1/3 to h = 2/3) between the asymptotic results 
(19) (22) and numerical results for the correlation function at a finite but large distance (£ — 30) is 
shown in Fig. 4. The numerical results are obtained by expressing the two-point correlator in the 
thermodynamic limit as the determinant of an £ x £ matrix (see section 3) and then evaluating the 
determinant for different times. As we are concerned with equal time correlators only we do not 
need to extract the two-point function from a cluster decomposition of the 4-point function [61]. 
The agreement is clearly excellent. The ratio between the exact numerics and the analytic result 
(19) in the space-time scaling limit is shown in the inset of Fig. 4 for two values of k = v max t/£. 
We see the ratio approaches a constant for large £. The corrections to this constant are seen to 
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be oscillating. The right panel of Fig. 4 shows results for the connected two-point function for the 
same quench (ho = 1/3 — > h = 2/3). We see that the numerical data are fit very well by both (19) 
and the form factor result (22). The connected correlator is exponentially small for t < t F and 
correlations start forming at tp. 

In order to elucidate the relaxational behaviour of the two point function it is useful to follow 
Refs [9, 22] and rewrite (19) in the form 

" r dkr 2t £ \ / 2t 

J v\V(kj~W)) H ^ 



p x p x p(t,t) 

- ' exp 



(23) 



Here we have defined mode-dependent correlation lengths and decay times r(fc) by 

C\k) = -ln|cosA fe | , T-\k) = -e' h (fc) In| cos A fe | . (24) 

We observe that these quantities are related by the velocity v k — e' h (k) of the momentum k mode 

v k T(k)=£(k), (25) 

which allows us to rewrite the theta- function in (23) as 9n(2vkt — £). The physical interpretation 
of (23) is then clear: a given mode contributes to the relaxational behaviour only if the distance 
£ lies within its forward "light cone" (the factor of two multiplying the velocity is explained in 
Refs [9, 22]). The form of the remaining factor then follows from the stationary behaviour: the 
time dependent piece compensates the factor (p FF ) , while the time- independent part is fixed 
by the t — > oo value of the correlator. As we already pointed out in our letter [71], this implies 
that the generalized Gibbs ensemble that characterizes the stationary state in fact determines the 
relaxational behaviour at late times as well. 



Approach to infinite times within the space-time scaling regime. For the quench within the 
ferromagnetic phase, the infinite time limit at fixed £ gives the same result as the infinite time 
limit within the space-time scaling regime, i.e. 



lim - In \ PF * F (l, t)\ = lim lim - In \ Pf * f (£, t)\ . 

t/i fixed 



(26) 



It is then useful to consider the approach to the stationary value within the space-time scaling 
regime result (19). In the limit t —> oo, we have 

dfc 



Pp x F (£,t — oo) oc exp 



In |cos Afcl 



(27) 



The corrections to (27) for large, finite times arise from the modes with e' h (k) ~ 0. For any h ^ 1, 
both modes with k — and k — tt contribute to this correction and at the same order since 
both ln|cosAfe| and the dispersion relation itself are quadratic at both points. A straightforward 
calculation gives for any h, h n < 1 



p x F x F (£,t^£) 



exp 



€/;^ln|cosA fc | 



oc 1 + 



(h-h Q ) 2 (l-2hh + hl) 



96tt(1 - h 2 f 



(^max^)' 



+ 



(28) 
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Figure 5. Numerical results for the connected two-point function for a quench from the 
ferromagnetic phase to a magnetic field very close to (left) and exactly at the critical point 
(right). The asymptotic behaviour agrees with the prediction (19) shown as continuous lines, but 
corrections are visible for smaller values of i. 



2.2.2. Quenches from the ferromagnetic phase to the quantum critical point. If we adjust the 
constant C§ F appropriately, equation (19) holds even for quenches to the quantum critical point. In 
Fig. 5 numerical results for the connected correlation function are compared with the asymptotic 
prediction for quenches very close to the critical point (left) and exactly to the critical point (right). 
In both cases, the asymptotic prediction is seen to become more accurate when £ is increased. 
However, it is clear from the figure that the asymptotic prediction works better for quenches exactly 
to the critical point than for quenches very close to it. This somewhat counterintuitive behaviour is 
readily understood as a property of subleading contributions to (19). We have already seen in the 
inset of Fig. 4 that there are subleading oscillating corrections. A more detailed analysis shows that 
they are power laws with an exponent that tends to zero upon approaching the critical point. This 
explains why the agreement of (19) with the numerical results is worse in Fig. 5 (left) than in Fig. 4. 
At the same time as the exponent of the power-law correction tends to zero upon quenching ever 
closer to the critical point, the oscillation frequency approaches zero as well. For quenches exactly 
to the critical point the leading oscillating corrections are therefore absent: they have morphed into 
a renormalization of the constant amplitude multiplying (19). This is the reason why (19) works 
better for quenches exactly to the critical point than for quenches close to it. 

For a quench to the critical point at h c = 1, the result (19) should be compared with conformal 
field theory (CFT) predictions of Refs [9, 22], which give the following result valid in the scaling 
limit of the Ising model 

p xx {r,t) J e -T*/8n, foYvt<r/2 



(l-h 2 )i X { e~" r /WvT " for vt > r/2 . (29) 

Here r is the physical distance, v the velocity characterizing the (strictly) linear dispersion relation 
e q = vq at the critical point, and To the so-called extrapolation time. Eqn (29) is valid for times 
and distances such that (x/v),t, (x/v) ~ t^> tq. The scaling limit of the result (19) is constructed 
in section 5, and considering a quench to the critical point we obtain from (264) 

lim 9 ^ \ oc exp [ -— \r9(2vt - r) + 2vt9(r - 2vt)} ) , (30) 
h-yi (1 - h 2 )i \ 2v I ij 
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Figure 6. Numerical data for a quench from the ferromagnetic to the paramagnetic phase. 
Left: The two-point function against the asymptotic prediction Eq. (32) (valid for t < tp) 
showing excellent agreement. Right: The connected correlation function for the same quench as 
on the left. For t < tp, p xx vanishes identically in the scaling regime. For t > tp, p xx (and so 
p xx ) approaches the asymptotic value Eq. (19), but subleading oscillating corrections are visible. 



where Ao = J(l — ho). Comparing (30) to (29) we see that the two expressions agree if we take the 
extrapolation time to be 

t = . (31) 

Outside the scaling limit the CFT expression (29) is not expected to provide a good approximation 
to the full result (19) because the nonlinearity of the dispersion relation becomes important. 



2.2.3. Quenches from the ferromagnetic to the paramagnetic phase. We have not been able to 
carry out a full analytical calculation of the time evolution of the two point function for a quench 
across the critical point. However, we conjecture! that (see Section 3.2.1) 

dk. 



P x F x P (e,t)~C£ P exp 



exp 



-2e' h (k)t\n\cosA k 



dk 



(£ - 2te' h (k)) In |cos A k \ H (2e' h (k)t - I) 



t > t F 



t < t f 



(32) 



+ cos(2e h (k )t + a) + . . . 

where a and k are the same as in (17) (see also Ref. [81]) and the constant factor Cp P is given in 
(18). This prediction is compared with the numerically calculated correlation function in Fig. 6 (left) 
and the agreement is clearly very good. For t < tp (32) is simply the square of the corresponding 
one-point function, which ensures that connected correlations vanish for t < tp in the space-time 
scaling regime. This is in agreement with numerical results for the connected two-point function as 
shown in the right hand panel of Fig. 6. As can be seen in Fig. 6 oscillations are present for t > tp 
as well, but they display a rather fast decay in time towards the determinant result (19). Like for 
the one-point function, the oscillating factor is a correction to the leading asymptotic behaviour in 



X Our conjecture is based on properties of the spectrum of the matrix T (50). 
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the space-time scaling limit, but it needs to be included to give a good description of the numerical 
data. Finally, we note that for h = 1 we have ei(fco) = and (32) reduces to (19). 

2.2.4- Quench within the paramagnetic phase. For quenches within the paramagnetic phase the 
form factor approach gives the following result (see Section 4.5) 

dk K{k) 



p x P x P (e,t) ~ p x P x P (e,oo) + (h 2 - i^VUVi r — sin(2te fe - kl) 

J-ir 1" £k 

- 2 f ^K 2 (p) {£ + 6 H (£ - 2te' p )[2te' p - £]) 
Jo 7r 



x exp 



+ 



(33) 



The regime of validity of (33) is sufficiently large values of I and t and "small" quenches in the sense 
discussed in the beginning of section 2. We have not attempted to calculate the infinite time limit 
pp x P (£, oo) within the framework of the form factor approach, because its exact large-£ asymptotics 
is known from the determinant approach to be [71, 81] 

p P x p(l,oc)~C£ P (t)e- e /t , (34) 

where Cp P (£) is determined in paper II [81] and 



£ = In (min[/i , hi}) — In 



hi 



h + ho 



2hh 



hi 



1 + hho 



l)(h 2 -l) 



(35) 



h + h 

As discussed in our previous letter [71], (34) is described by a general Gibbs ensemble. Based on 
the form factor result (33) one may speculate that the full answer may have the structure 

dk K(k) 

£k 



P xx (£,t) 



C^(£) + (h 2 



l^VU 2 ^ f - K( ^ S m(2te k -k£) + ... 

J-tt K ?k 



x exp 



dp 



In 



l + K 2 {p) 
1-K 2 (p) 



{£ + e H {£-2te' p )[2te' p -£]) 



(36) 



In Fig. 7 we compare the analytic result (33) to numerical results obtained for two different 
quenches within the paramagnetic phase. The agreement is seen to be excellent. 

For strong quenches the form factor result is not expected to be quantitatively accurate. This 
can be seen in Fig. 8. In all cases, the two-point function is seen to display slowly decaying oscillatory 
behaviour on the time scales shown. At sufficiently large t the decay is proportional to t~ 3 l 2 . This 
is in marked contrast to quenches within the ordered phase. The origin of this difference lies in the 
nature of the relaxational processes that drive the time evolution. The oscillatory behaviour seen 
in the paramagnetic phase arises from processes involving the annihilation of spin-flip excitations, 
while the smooth exponential behaviour seen in the ferromagnetic phase is related to the ballistic 
motion of domain wall excitations. 

The structure of (33) implies the existence of a late time crossover scale t* , at which the second 
contribution becomes smaller than pp x p (£,oo). Using that pp x P (£, oo) oc e^ e ^ and that the second 
contribution decays like ( Jt)~ z / 2 at late times we may estimate t* as 

Jt * „ e vm) . (37 ) 

For the cases shown in Fig. 7 this gives Jt* ~ 4114 and Jt* ~ 10 20 respectively (Jt F = 7.5 in both 
cases). This means that in both cases the stationary behaviour characterized by the generalized 
Gibbs ensemble is revealed only at very late times. 

Finally, we note that so far we have not been able to analyze the time evolution of order 
parameter correlators for quenches within the paramagnetic phase by means of the determinant 
approach. 
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Figure 7. Time evolution of order parameter correlators for quenches within the paramagnetic 
phase. Left: pp x p (l = 30, t) for a quench from h$ = | to h = 2. Right: pp x p (£ = 30, t) for a quench 
from ho = oo to h = 5. The two-point function is seen to exhibit oscillatory power-law decay 
at late times. The form factor result (solid lines) is seen to give a very good description of the 
numerical data (points). The short time regime is not shown as the correlators are exponentially 
small by virtue of the horizion effect. 




Figure 8. Time evolution of p P x P (£ = 60, t) for a quench from ho = 5 to h = 5/4. As the quench 
is no longer small, deviations from the form factor result become more pronounced. 



2.2.5. Quenches from the paramagnetic to the ferromagnetic phase. Here we have not been able to 
obtain analytic results by either the determinant or the form factor approach. However, observations 
within the framework of the determinant approach (see Section 3.2) suggest that for late times t > tp 
the leading asymptotic behaviour of the two-point function should be given by 

— \n\cosA k \0 H (2e' h (k)t-i) 



^(£,i)~Cj5 F (£)exp 



x exp 



2t 



71 dk 



e' h {k) ln|cosA fc |0 H (/-24(A;)*) 



t > t F . 



(38) 
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Figure 9. Order parameter two point function for a quench from the paramagnetic to the 
ferromagnetic phase (inset) and its ratio with the asymptotic prediction ppp given in eqn (38) 
(main plot). We consider a quench from ho = (2 + 3\/2)/2 to h = 1/2. For any t > tp, the ratio 
approaches 1 with increasing I. 



Here the factor Cp F (£) is given by the t — oo results of [81] 

V(i-ft 2 )(ft§-i) 



CpfW = 



ho — h 



cos ( £ arctan 



1 + h h 



(39) 



We have tested this conjecture by comparing it to the numerical results and found it to hold. An 
example is shown in Fig. 9, where we plot the ratio of the numerically calculated correlation function 
and the analytic expression (19) for several values of the distance i. The ratio clearly approaches 1 
at late times. We note that the values of I have been chosen in a way such that for the particular 
quench considered (i.e. h = 3.12 to h — 0.5) oscillations in I are suppressed in the t — > oo limit. 
Our analytic methods do not currently provide an understanding of the quench dynamics for times 
shorter than the Fermi time t < tp. 



2.3. Quenches to h — or h = oo. 

These quenches are special, because the post-quench Hamiltonian is classical in both cases. As 
a result the dynamics is anomalous. For example, for quenches to h = correlation functions 
involving only erj operators are time independent, because the final Hamiltonian H(h) commutes 
with all o~j. Correlation functions involving other operators do depend on time, but generally 
exhibit persistent oscillations. For instance, we have that 

e iff(fc=o) to .. e -ur(fc=o)t = cos 2( 2 jt) CT | + I sin ( 4 Jt) + of +1 )o-V - sin 2 (2Ji)aJ_ 1 a|aJ +1 , (40) 
and as a result (t)\a z ^ (t)) does not approach a stationary value at late times. 

3. Determinant approach: Analytic derivation of the asymptotic two-point function 
for a quench within the ordered phase. 

In this section we present the analytic derivation of Eq. (19) for the asymptotic behaviour in the 
space-time scaling regime (t, £ 3> 1, but t/£ arbitrary) of the two-point correlation function p xx (£, t). 



Quantum Quench in the Transverse Field Ising chain I 



16 



Within this approach, it is convenient to replace the fermions Cj in Appendix A with the Majorana 
fermions 



a*=c\+Ci a v -=i(c\ 



"j — "j ' ^3 ">j — "Vj ~~ yj) > (41) 

which satisfy the algebra {af, a*} = 2<5;„, {af , a^} = 2Si n , {af , a^} = 0. In terms of these Majorana 
fermions, the operator <rj has the nonlocal representation 

t-i 

< IT (42) 

The two-point function of <r x is then the expectation value of a string of Majorana fermions 



p™(e,t) = (H(-iay(t)a] +1 (t))). 
The Ising Hamiltonian (1) can be written as (see Appendix A) 



H 



1 



IT 



H 



1 



i=i 



(43) 



(44) 



where H and _ff e are quadratic Hamiltonians of the Jordan- Wigner fermions in the sectors with 
odd and even fermion number respectively {H Q / e both commute with Y\^ =1 erf). For finite chains, 
the ground state of the Ising Hamiltonian is an eigenstate of H e . However, in the thermodynamic 
limit, the Z2 symmetry of H is spontaneously broken in the ordered phase ho < 1, where the ground 
state is a linear superposition of the ground states of H e and H a 

|B(0))ns + \B(0)) R 



l*(0)> 



(45) 



Thus, for a quench starting in the ordered phase free fermion techniques cannot be straightforwardly 
applied to the calculation of correlation functions involving generic operators. However, here we 
are interested only in the expectation value of even operators O e characterized by 



















.1=1 . 



= o e . 



(46) 



Hence we have 



lim (*(t)|O e |*(t)> = lim 



NS 



(B(t)\d e \B(t)) NS + R (B(t)\O e \B(t)) K 



L— >oo 



L— >oo 

lim m (B(t)\O e \B(t)) NS 

L— J- 00 



(47) 



Crucially, the expectation value Ns(B(t)\O e \B(t))^s can be evaluated using Wick's theorem. In 
contrast to even operators, expectation values of odd operators 



















.1=1 



-On 



(48) 



are significantly more difficult to determine [62]. The particular case of interest, p xx (£,t), involves 
an even operator (cf. Eq. (43)) and by straightforward application of Wick's theorem one obtains 
a representation as the Pfaffian of a It x 2t antisymmetric matrix [62] 

p xx (t,t)=pf(T), (49) 
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where pf denotes the Pfaffian, T is given by 

r r_i ••• ri_£ 
ri r ■ 



-fl 9i 
~9-l fl 



(50) 



r^-i r 

The elements of this matrix are the fcrmionic correlations!) 



9n 



l f a l+n-l) 



in + iS n o = i (afaf+n) 



i {a v i +n a v i ) 



(52) 



The matrix f is a block Toeplitz matrix because its constituent 2x2 blocks depend only on the 
difference between row and column indices. It is customary to introduce the (block) symbol of the 
matrix T via Fourier transform as follows 



e uk T(k), with r(fc) 



-fl 9i \ = f v dk 
-9-i fl ) J_ v 2tt 

where the functions f(k) and g(k) are 

f(k) =ismA k s'm(2e h (k)t), 

g (k) = -e iek ~ ik cosA fe - isinA fe cos(2£, l (fc)f) . 



-/(*0 g(k) 

-9{-k) f(k) 



(53) 



(54) 



The spectrum of block Toeplitz matrix f is the same of an £ x I Tocplitz+Hankel matrix (Hankel 
matrices have elements which depend only on the sum of row and column indices, instead of the 
difference as in Toeplitz matrices) . Indeed, being T a real antisymmetric matrix, its eigenvalues are 
complex conjugate pairs ±i\j. Given an eigenvector V of T with eigenvalue iX, we can define the 
vector w with I components w; L = Vii-\ (since r_ ; = a y YiG y , taking the even components would 
change the sign of A). Thus, the i vectors Wi are solutions of the eigenvalue problem 

{iT ± H)w, = TiXjWj, { 2: = J l ~ n T » (55) 

^ 11 ln — yl+n— L — 1 

In particular, this means that 

|pf(f)| = |det[iT±if]|. (56) 

The sign of the Pfaffian can be fixed by observing that, for any given £, both the determinant of the 
Tocplitz+Hankel matrix and the Pfaffian of the block Toeplitz matrix are polynomials in /; and gi 
of the same degree. Thus, the sign is independent of the actual values of the matrix elements and 
can be determined by considering the simplest case T = (which occurs at the initial time and at 
asymptotically late times after the quench), obtaining 



(-1)*-^ det(H + iT) 



(-1)— s— detW, 



(57) 



§ The identity in the second line of Eq. (52) is a consequence of reflection symmetry, which indeed implies that 
correlations should be invariant under the transformation 

af-xQl^H+i-! < ^ -'(FKH+1-! ■ ( 51 ) 
I I 

For example, the Dzyaloshinskii-Moriya interaction breaks this symmetry. 
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where W = H + iT. The matrix f is real and antisymmetric, so its eigenvalues are complex 
conjugate pairs of the form with j = 1. . .£. By construction Xj € [—1, 1] (see e.g. [82]) and 

they are also the eigenvalues of the matrix W = H + iT. Thus we have 

e 

{p xx f = det f = dot W 2 = Yl X} • (58) 

j=i 

The evaluation of the correlation function p xx (£, t) for large I is then equivalent to the 
asymptotic evaluation of the determinant of a block Tocplitz matrix or the sum of a Toeplitz 
and a Hankel matrices. Such matrices have been the subject of intense study by mathematicians 
and physicists for more than a century and standard, rigorous techniques for calculating their 
determinants such as Szego's lemma and the Fisher-Hartwig conjecture are available, see e.g. 
[83, 84]. However, these methods have been specifically designed for the evaluation of determinants 
of matrices whose elements do not depend explicitly on the matrix size. This is in contract to 
our case, where we are interested in the scaling limit t , t — > oo with finite ratio £/t. Under 
these circumstances, each clement of the matrix T in the scaling limit depends on a parameter 
(namely t) which is proportional to the matrix dimension 21. This precludes the application of the 
aforementioned techniques for the asymptotic evaluation of these determinants. The two exceptions 
are the t = case, where we recover the known equilibrium results, and the limit t = oo. In order 
to deal with large values of t we developed a completely novel approach, which follows the one we 
proposed in Ref. [66] for the entanglement entropy and it is based on a multi-dimensional stationary 
phase approximation. 

In the following we derive a rather general result valid for any block 2x2 block Toeplitz matrix 
T with a symbol t(k) that can be cast in the form 

t(k) = n x {k)a x k) + n ± (k) ■ 3™ e 2le{k)ta * k) , n ± {k) ■ x = . (59) 

Here the time t is the only parameter proportional to the matrix size 21, n x , n± are fixed but 
otherwise arbitrary and denotes a local rotation of the Pauli matrices 

a (k) _ e Lw(ky<? aae -iw(k)-a (g ) 

Our particular case of interest (53) corresponds to having n x + \n±\ 2 = 1. 

We first consider TrT for positive integer n. We note that the traces of odd powers of 
r vanish, because T is a real antisymmetric matrix (r is diagonalizable and for each non-zero 
eigenvalue there is another one with opposite sign). Our main result is that 



Tr[r 2 »] _ r ol ^,„ „t 

t/e fixed 



t _lim / ^max(l-2| £ '(fc o )|-,0)(n ;c (fc o ) 2 + |n ± (fc o )| ; 



+ f ^min(2| £ '(fc )|^l)n ;c (fco) 2 ". (61) 

Given the result (61) it is possible to infer the asymptotic behaviour of more complicated quantities 
such as 

Tr^(r 2 ) = J2 Tr[r 2n ] , (62) 

n 

where F(z) is an analytic function with power series expansion J-{z) = J2 n FnZ n around z = 0. 
Using Eq. (61) and interchanging the order of sum and integration, we have 

Tr[j-(r 2 )] _ r dfcc, / o . ..t 



t,i- 

t/l fixed 



M h = / ^max(l-2|e'(fc o )l^0)^(fc o ) 2 + l^(fco)| 2 ) 
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+ f ^min(2| £ '(fc )|*,l)j-(n ;c (fc ) 2 ). (63) 

Clearly this approach is fully justified only as long as all eigenvalues of T fall within the radius of 
convergence of the expansion of the function F{z) around z = 0. For example, in the case of the 
entanglement entropy, we have T[z) = — ^ip 2 In [66] and since the eigenvalues of T are of the 
form with Xj G [—1,1] (63) holds (further generalizations to the entanglement of two blocks 

have also been considered [85]). 
In the case at hand we have 

(p xx ) 2 = detr = exp[Trlnr] , i.e. ln^) 2 = X - TrlnT 2 . (64) 

In order to use (62), (63) we are therefore led to consider the function T{z) = |lnz. The latter 
has a branch point at z — and the previous approach appears not to be applicable. In order to 
circumvent this problem we employ a power series expansion of the logarithm around z = 1 

ln(p xx ) 2 = - Tr In T 2 = - Tr ln[l + (r 2 - 1)] = - V ^ '- Tr[(r 2 - l) m ] . (65) 

A A A Tfi 

m— 1 

Application of (63) to the function T{z) — (z — l) m then results in 

Jm Tr[(r2 2 ; = r < ^max(l-2| e '(fc o )l^0)K(M 2 + l^(fco)| 2 -l)" 1 

t/t fixed -'-■n" 

+ f ^min(2| £ '(fc )l*,l)(^(M 2 -ir- (66) 
Finally we can carry out the sum over m 

Jm * r2] = r^max(l-2| £ '(fco)l',0)ln(n x (fc o ) 2 + l^(M| 2 ) 

t/t fixed •> 

+ f ^min(2| £ '(fc )|^l)ln(n :c (A :o ) 2 ), (67) 

which is exactly the same result we would have obtained applying (63) directly to the function 
F{z) = \h\z. The reason for this lies in the simple algebraic dependence of (61) on n. The 
specialization of (67) to the longitudinal correlator p xx is now straightforward. The symbol f (A:) 
of the block Toeplitz matrix T can be cast in the form (59) with 

w(k) = [6{k) - k]x , n 2 (fc) = cos 2 A fe , |n_L(fc)| 2 = sin 2 A fe . (68) 

Eq. (67) then gives the asymptotic behaviour 

t lim jln\p xx (£,t)\= ^ ^min(2| £ '(fc )l^l)m|cosA fc0 |. ( 69 ) 

t/t fixed J— IT 

This is a rewriting of (19) and represents one of our main results. 

In the reasoning leading up to (67) we have assumed that the eigenvalues of T 2 — 1 lie within 
the unit circle so that we can expand the logarithm in a power series. This is equivalent to the 
requirement that the eigenvalues Xj of T are such that < A 2 < 1. More precisely, if there exists a 
x > such that for all i the eigenvalues Xj of T fulfil 

z < A 2 < 1 , (70) 
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then the results (67) and (69) hold. In some cases of interest we find that even though Xj are 
different from zero for any finite £, one or several eigenvalues approach zero in the limit £ — > oo. In 
these circumstances the steps leading up to (67) and (69) are no longer justified. However, even 
then (69) can provide useful information about the asymptotics of p xx (£,t). On general grounds || 
we expect the latter to be of the form 

lim \n\p xx (£,t)\ = -£p(t/£)-a\n(£) + c(t/£) + ... , (71) 



t/£ fixed 



and the question is under what circumstances p{t/£) is given by (69). Under the assumption that 
the general structure (71) holds, this depends on the number No of eigenvalues Xj that approach 
zero for £ — »■ oo, as well as on how quickly they tend to zero. If N is finite and the corresponding 
eigenvalues approach zero sufficiently slowly with £, (69) will still be applicable. As an example let 
us consider there case where only a single pair ±iXi of eigenvalues approaches zero for £ — > oo in 
such a way that for large £ we have |Ai| ~ £~@. Since p xx (£,t) — . X 2 , the eigenvalue pair will 
contribute to subleading logarithmic term in (71), but will leave the function p{t/£) unchanged. On 
the other hand, if our pair of eigenvalues were to approach zero exponentially fast |Ai| <~ e - ^*/^, 
then it would contribute additively to the function p{t/£) and (69) would cease to hold. 

As far as our quench problem is concerned we do not have a criterion that would establish a 
priori whether eigenvalues exponentially close to zero will be present. We find that for quenches 
within the ferromagnetic phase they arc always absent. In the other cases, the picture emerging from 
numerical studies of the spectrum of f suggests that only for quenches starting in the disordered 
phase eigenvalues approach zero exponentially fast, and then only a single pair ±iAo does so. 
Interestingly, we find that the contribution of all other eigenvalues is again captured by (67). At 
present we are not able to determine A analytically. 

3.1. Proof of the formula for the trace of integer powers ofT. 

The first step of our calculation is to derive an appropriate integral representation for TrT 2n by 
trading matrix multiplications for additional integrations. The basic idea is easily explained for a 
product of two r matrices. Each block element is given by 

Tin = f ^ l - n H{k) , (72) 

where the 2x2 block t(k) is of the form (59). Hence the block matrix elements of T 2 can be written 

as 

(r 2 ), m = J2 T ^T mn = I ^^e^-^Hih) ■ i(k 2 ) J2 e im(fel - fc2) • (73) 

m=l L ">"J m=l 

ow be replaced by an integral using 

i)fe/2 V e" nk - - C (\f k r tlik/2 (74) 
hx "2 7_ 1 ^2sin(fc/2) e ' ^ 

The generalization to T 2n is straightforward, and replacing the trace in an analogous way we obtain 
the following integral representation 

Tr[r 2 «] = (^) 2 "| ^ J d^CtttyeH^W^-WFUk}). ( 75 ) 

[ — 7T,7r] n [-14] 2 " 

|| This form has been observed is various examples, see e.g. [22, 34, 66, 61, 77, 36, 56], even with oscillating factors. 



The sum over m can now be replaced by an integral using 

-i(e+ 



e 
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Here the functions appearing under the integrals are 

2n-l 



^(«) = n 



/i 2sin[(fc,-fc,_ 1 )/2] ' 



2n-l 
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(76) 
(77) 



(78) 



= Tr[ [] M*i)4 fc) + fix(fei) • iW"*^ 

i=0 

the trace is over the remaining 2x2 blocks and £j = e(fcj). Performing the change of variables 

Co 

C* = Ci+i — «e[l,n-l] 
Eq. (75) can be rewritten in the form 

Tr[r 2 "] = J ^ J d 2 "C C({fc}) e-^ZT^to-ko)/* F({fc}) . (79) 

[-7T,7r] 2 " -R C 

Here the domain of integration i?j is determined by the conditions 

fc-i 

i? c : -1 < ^ Cj < 1 Vfc e [!) 2n l • ( 80 ) 
j=o 

As the integrand in Eq. (79) is independent of Co, we can carry out the Co integration, which gives 



£\ 2n 



Tr[r 2 "] = (£) 



d 2n fc 



(27T) 



2// 



d 2 " _1 C M({C}) C({fc}) e-^^i'Ofe-feo)^ (si) 

[-ir,7r] 2 " 

Here the function /u({C}) is the measure of the domain of Co under the constraints (80) 



M({C}) = max 



0, min 

je{0,2n-l} 



fc=l 



mm 

ie{o,2«-i} 



1 



fc=i 



(82) 



It can be shown that ^({C}) is symmetric with respect to an arbitrary permutation of the variables 
Ci- Since we are interested in the behaviour for t ^> 1 and the phase in the integral is proportional 
to the large parameter £, the asymptotic behaviour can be obtained using a multi-dimensional 
stationary phase approximation. The main idea behind this method is that the leading contribution 
to the integral arises from the neighborhoods of the points in which the phase is stationary. As the 
symbol is independent of the integration variables Ci, the stationarity of these variables implies 

kj^ko, j = l,...,2n-l. (83) 

We may replace any kj with fco everywhere except in rapidly oscillating terms such as the e 2teit 
factors in the symbol. We call this the localization rule. This rule allows us in particular to drop 
the factor C({fc}), cf. (76), as it is equal to one at the stationary points C({fc}) ~ 1 



Tr[r 2 



£,2r, 







d 2 "fc 

(2tt) 2 " 



L C 



2 F({k}). 



(84) 



The localization rule furthermore allows us to substitute the integration variables in n x and ri± 
with k 0} and to remove the local rotation in /({&}), i.e. under the integral we may replace 

2n-l 



F({k}) -+ G({k}) = TrJ[ [n x (k Q )a x + n ± (k ) ■ Se 



2isitcr x 



(85) 



i=0 
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As is shown in Appendix C, the product under the trace can be rewritten as 



2n 

J2n 2 x n - p \n ± \P{ia z )P ^ (-ijEJUi im e «to. E^=i(-i) p - m ^ m _ 1 . ( 86 ) 

p=0 l<ji<j2<-..jp<2n 

Each term in the sums over ji . . . j p gives the same contribution to the integral (84) as can be shown 
by changing integration variables 

^1 = kji > ^2 = kj2 1 ' ' ' i — ^pi k[ — ki for all Z ^ jji , . . . , jp} , 

Ci - Cn, C 2 = Cj 2 , ' ■ • ,C jp = C P , Ci = Ci for all I $ {j u . . . ,j p }, (87) 

and then invoking the invariance of ^({C}) in Eq. (82) under any permutation of its variables. 

When evaluating the trace in (86) we may therefore replace j m — > m. We call this the contraction 

rule. Noting that only the terms with even p have a non- vanishing trace (because Tr e la<Tx = 0) 

and then applying the contraction rule allows us to replace 

n 2p 

G({k})^J2( I )nT- 2p n 2 f2cos(2tJ2(-ire m ^). (88) 

p— m— 1 

Here the binomial takes into account the number of terms giving identical contributions to the 
integral and we have used that Tr e ta<Jx = 2 cos a. Inserting this expression into (84) we arrive at 

x nl^nfe-^T^ 1 O^+^SVD^ , (89) 

where we replaced the cosine with a complex exponential using the symmetry of the integral under 
the change of variables k — »■ —k and ( — > — (. 

We are now in a position to employ a stationary phase approximation in order to extract an 
exact asymptotic result in the limit of large I. As the phase of the integral is stationary along a 
one-dimensional smooth curve, application of the stationary phase method is not a simple matter. 
For a two-dimensional integral whose phase is stationary on a one-dimensional variety, the solution 
can be found in Ref. [86] . Higher-dimensional integrals with a curve of stationary points are only 
partially treated in Ref. [87], where it is demonstrated that it is possible to isolate the integration 
in ko and perform a standard multi-dimensional stationary phase approximation for the remaining 
integrals. To apply this idea to our case, we rewrite Eq. (89) as 



Trlr-H^) 2 "" 1 ^ " ) f_^n x {k r- 2 W{k Q rK,{k a ), 



(90) 



1=0 

where the functions A„ ; ; are given by 



K-M = J |^ / d-^c M{C}) e-^r^c^wf- ^ . (91) 

[-tt.tt] 2 ™- 1 

We now use the standard multi-dimensional phase approximation to evaluate A n -i(ko). 

The general result for the large-i? asymptotic behaviour of a multi-dimensional integral over 
rapidly oscillating functions, whose phase is stationary at an isolated point x detached from the 
boundary is [86] 

d N xp(x) = ( 2 J) Ar/2 p(f )|detA|- 1 /2expU(f ) + ^l , (92) 



i ) 
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where A^ — d Xi d Xj q\% is the Hessian matrix of f(x) evaluated in xq, and a a is the signature of 
the matrix A, i.e. the difference between the number of the positive and negative eigenvalues. The 
stationarity conditions for the phases of the integrands in A„ ; ;(fco) are 

kj=k , j = l,...,2n- 1, 

C j= 4|(-1)%;., j = l,..., 21-1, (93) 
=0, j = 2l,...,2n-l. 

Using an ordering of integration variables where the £'s are placed before the fc's, the Hessian is of 
the form 



2 V 1 M , 

Hence the eigenvalues cl± of A are related to the eigenvalues fii of the matrix M by = 
{Hi ± y / /xf~+4)/4. The signature of A is = and its determinant is det A = -4 1 - 2 ". At 
the stationary point the phase in the integral vanishes, because there is an even number of e's with 
alternating signs. Finally, the value of /i({C}) a t the stationary point is found to be 

M ({C}) = | |max[0,£-2| £ '(fc )|i] j^O, (Q5) 

Putting everything together, the stationary phase approximation gives the following result for the 
functions A n .;(fco) 

A niI (fco) = 0) M{C»- (96) 
Inserting (96) into (90) then gives 

Tr[r 2 «] =iJ2{ ^ ) /^-,(M 2 "- 2 V(M 2 V({C}), (97) 

and substituting the value (95) for /x({C}) we arrive at 

Tr[r 2 "]=2 y" ^max(£-2|e'(fc )|i,0)(n ;c (fc ) 2 + |n ± (fc )| 2 )" 

+2 ^ ^°^-max(0,^-2| £ '(A : o)|t))n :c (A : o) 2n . (98) 
This is equivalent to (61). A less general version of this result was previously presented in Ref. [66]. 



3.2. On the applicability of stationary phase method. 

Eq. (69) is our main result obtained with the determinant approach. It is based on a multi- 
dimensional stationary phase approximation for Tr[r 2 ™]. The result (98) is then used in combination 
with a series expansion in order to obtain the asymptotic behaviour of det T, which in turn gives 
the square of the longitudinal correlation function p xx . As we have already discussed, this series 
expansion is possible only if there is no eigenvalue of the matrix T that approaches zero in the 
large I limit. While this restriction appears to be quite simple, we do not have an analytic method 
that allows us to predict for what kind of quenches zero eigenvalues exist in the t — > oo limit. To 
address this question we therefore have carried out numerical studies of the spectrum of W (which 
is related to the spectrum of f by (58)) for different quenches. We have to distinguish between four 
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Figure 10. Typical absolute values of the low lying eigenvalues of the matrix W in the four 
kind of quenches. In (a) and (b) I = 30, while in (c) and (d) I = 60. (a) For a quench within 
the ordered phase, all eigenvalues are well separated from zero. Before the Fermi time (up to 
finite distance corrections to tp ) the eigenvalues are almost degenerate in pairs, (b) For a quench 
within the disordered phase, a single eigenvalue is very close to the real axis (exponentially close 
for t < tp and as a power law for t > tp). (c) Quench from the disordered to the ordered phase. 
For t < tp a single eigenvalue is exponentially close to the real axis. For t > tp the distance of the 
smaller eigenvalue from the real axis scales like a power law. (d) Quench from the ordered to the 
disordered phase. All eigenvalues are doubly degenerate. All the smallest eigenvalues approach 
zero as a power in I for any time. For t < tp the smallest eigenvalue crosses periodically zero. 



cases: (1) FM — > FM; (2) PM — > PM; (3) PM — ► FM and (4) FM — > PM, where FM and PM 
denote the ferromagnetic and paramagnetic phases respectively. In Fig. 10 we report the (absolute 
values of) smallest eigenvalues of the matrix W for particular examples of the four cases. Extensive 
numerical studies suggest that the spectra are similar for all quenches of a given type, i.e. (l)-(4). 
The qualitative features emerging from Fig. 10 can be summarized as 

(1) FM — >■ FM: all eigenvalues are well separated from zero; 

(2) PM — > PM: a single eigenvalue of W is close to zero. For t < tp it approaches zero 
exponentially fast in £, while for t > tp it tends to zero only as a power-law; 

(3) PM — > FM: for both t < tp and t > tp several eigenvalues scale to zero like power laws in £. 
For t < tp the smallest eigenvalue approaches zero exponentially fast in t. 




Figure 11. Ratio between the numerically evaluated two point function p xx (£, t) and the smallest 
eigenvalue of the matrix W. Left: Numerical data for a quench within the paramagnetic phase 
from ho = 3 to h = 3/2. The data are in perfect agreement with the prediction (69) and the 
leading correction is a stationary multiplicative factor. Right: The same plot for a quench from 
the paramagnetic ho = (2 + 3\/2)/2 to the ferromagnetic h = 1/2 phase. 



PM: several eigenvalues tend to zero in a power-law fashion in i for both t > tp and 
For t <tp the smallest eigenvalue crosses zero periodically. 

The numerical analysis suggests that (69) can be applied for all quenches of type (1) for large t and 
all corrections to (69) are small. In case (4) our result (69) still gives the dominant behaviour, but 
there are additional important oscillating power-law corrections to p xx (£,t) that we have not been 
able to determine analytically. For quenches originating in the paramagnetic phase (69) gives the 
dominant contribution only if t > tp, but there are important power-law corrections to p xx (£,t) 
beyond the accuracy of our analysis. 

These properties of the spectrum suggest that, for quenches starting from the disordered phase, 
the product of all eigenvalues except the only one which is close to zero (i.e. p xx /Aq) should be 
given by Eq. (69). Fig. 11 shows that indeed this is the case for quenches to both disordered and 
ordered phases for any t. Thus, if we were able to obtain a prediction for the smallest eigenvalue, 
we would completely characterize these quenches as well. 

Finally, we note that the gross features of the spectrum of W (and hence T) appear to be 
related to the winding number (around zero) of the (block) symbol (53). We observe that the 
winding number of g(k) (which is only a part of T(fc)) depends on the quench parameters in the 
following way: 

(1) If \h\, \ho < 1, the winding number is zero; 

(2) if \ho\, \h\ > 1, the winding number is equal to 1; 

(3) If \h \ < 1 and \h\ > 1, the winding number oscillates with t between and 1. In particular it 
remains for any t > tp. 

(4) If \ho\ > 1 and \h\ < 1, the winding number oscillates with t between —1 and 1. 

We see that the winding number vanishes when there is no eigenvalue exponentially close to zero, 
but the reverse does not hold. For (block) Toeplitz matrices with £ independent elements the 
generalization of the Szego lemma to symbols with non vanishing winding number is complicated 
and not known in general (see e.g. [83, 84]). 



(4) FM — > 
t < t F . 



Quantum Quench in the Transverse Field Ising chain I 



26 



3.2.1. A closer look to quenches from the ordered to the disordered phases. For quenches from the 
ordered to the disordered phase, the last panel of Fig. 10 shows that all eigenvalues of W move 
coherently in time. In particular, the eigenvalues closest to zero oscillate with the same frequency 
about zero until the Fermi time t F . After tp they do not cross zero anymore and the smallest one 
remains separated from zero as a power law in L These qualitative features of the spectrum would 
suggest that for t < tp p xx (t,£) should oscillate in time with the same frequency of the smallest 
eigenvalues. While it is not easy to put this on a solid ground, we have been able to identify it as 
the frequency with which the winding number of g{k) changes. We indeed find that the frequency 
is determined by the mode ko such that cosAfe = 0. Thus for t < tp we heuristically predict 

PIT (^fe 

(af(t)af +1 (t)) ~C£p[l + cos(2£ feo £ + Q!)]exp[y ^2|e' fe |iln | cos A fe |J t < tp , (99) 

where a is an undetermined phase shift. The comparison of this prediction with the direct 
computation of the determinant has already been shown in Fig. 6. For times larger than tp, 
Eq. (69) matches well the evolution of the two-point function, up to oscillatory decaying corrections 
(see Fig. 6) which is consistent with the presence of one eigenvalue exponentially close to zero. 

4. Form-factor approach 

While the procedure set out in Section 3 can be generalized to any model with a free-fermion 
representation [88], it cannot be applied to quenches in interacting integrable models. In order to 
overcome these limitations we have generalized the form-factor approach to correlation functions 
in integrable quantum field theories [72, 73, 74, 75, 76] to quantum quenches. A characteristic 
feature of integrable models is the existence of a basis of scattering states of elementary excitations, 
which are simultaneous exact eigenstates of the Hamiltonian and the momentum operator. These 
states can be characterized by n momenta k\ , . . . , k n and in general also other quantum numbers 
a 1 ,. . . ,a n 



H\ki, . . . , k n ) aii ,,, t a n — |^ ^ ' £q ; - jkj) \ki, . . . , k n ) aii 

3 = 1 

71 

P|A;i,...,fcn)oi,...,o n = \ [ ^k^\k 1 ,. . . ,k n ) ai ...., an . (100) 

Ground state correlators of local operators O can be expressed in a Lehmann representation 

oo 

(O|0(t,aO0t(o,O)|O>= ^ l(0|O(0, 0)|fc 1; . . . , fc„) ai ,..., a J 2 

n— fci,...,fc n 
ai,...,a n 

X e -*5:"=l e aj(*i>*-*J* } (101) 

where |0) denotes the ground state and in our notations corresponds to particle number zero. 
Using the known expressions for the form factors (0|O(0, 0)|fc 1; . . . , fc 7 i) ail ... j0n the spectral sum can 
be calculated to very high accuracy (for large t and x) by taking into account only terms involving 
a small number n of particles. Recently the form factor approach has been generalized to low 
temperature correlation functions [89, 90] (see also [91, 92] for related work) 

1 Tr[e-^] (1 ° 2) 
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The numerator in (102) can now be expanded in a Lehmann representation as follows 



E 



n,m— ki,...,k n 



E 



6i(Pr7 



, Pl \O(0,0)\k u ...,k n ) au 



x e 



(103) 
[72, 93], one 



Using exact results for the form factors b m ,...M (Pmi ■ • • ,Pi\0{®, 0)|fc 1; . . . , fc„) aii .. 
needs to sum an infinite number of terms in the Lehmann expansion in order to obtain the correlation 
function at late times and large distances. Such a resummation is possible at low temperatures 
T <C A, where A is the spectral gap, because the density of excitations in the state of thermal 
equilibrium constitutes a natural small parameter in this case [89]. The situation after a quantum 
quench bears many similarities to the finite temperature case. When dealing with a "small" quench 
in a gapped theory there exists a regime in which the density of excitations (of the post-quench 
Hamiltonian) in the initial state is small. It is then natural to use this small parameter in order to 
carry out a low-density expansion for the observables of interest. Let us consider a quantum quench 
H a — > H, where H describes an integrable scattering theory: we prepare the system in the ground 
state |\&o) 01 Ho and then consider time evolution by H. We are interested in observables such as 

(*o(t)\0(x)Ot(p)\* (t)) 



<*o|*o) 

Any translationally invariant initial state can be expressed in terms of the eigcnstates of H as 



(104) 



l*o(*)> = £ E e-^^ k ^F n ({k 3 ,a 3 })\k u 

n—0 fei fen 

ai ,...,a n 

F n ({kj, aj }) = 0n ,..., ai (k n , . . . , fci|tf (0)). 
A particular class of initial states is given by 



■ • i &n)ai,.. 



i*o(o)HE^y E 



n=0 pi, 



U k " j " j m 



Pi, Pi, 



-Pn i Pn) ai ,&i ,. . .,a n ,6 n 



(105) 



(106) 



where summation over the indices aj, bk is implied. As is shown in Appendix B quantum quenches 
of the transverse field in the TFIM automatically lead to initial states of the form (106). For general 
integrable quantum field theories (106) are known to describe boundary-states [94] and correspond 
to situations where the initial state can be viewed as a boundary condition (in Euclidean space) 
that is compatible with quantum integrability. Given an initial state of the form (106), the basic 
idea is to employ a Lehmann representation in terms of the exact eigenstates of the post-quench 
Hamiltonian H. For a two-point function we have 



(Mt)\o(x)oHo)\Mt)) = 



X - x - I 



E 



n!m! 

m=0n=0 Pi,-.-,Pn 



IK" ' ■:/',: 



3 = 1 



Y[K c ' d '(kt 



><E E 

s=0 gi,...,q 



1 

^! c 

■ Ss,- 



% 4 m ,..., Cl 4i{ k rn, -k m , ■ ■ ■ ,k\ - fci \0(t, x) \q x , . . . ,q 8 )f u ...,f a 

jMs, ■ ■ -,<ll\0 ] {t,Q)\ - Pi, Pi,.. .,-p„,P„> 0ll 6 1 ,..,o n ,6„- 



(107) 
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The form factors entering this expression are known, both for the TFIM and a variety of intcgrable 
quantum field theories [72, 93]. The normalization itself has the following Lehmann representation 

<*o(t)i*o(*)>=£w- j2 n^'fou^'to)]* 

n=0 ^ Pi.---.Pn ,' = 1 

fci,...,fc n 

X c„,d„,...,ci,<ii(fcn,-*n,---,fel -Pl,Pl,---,-Pn,Pn>oi,6i,...,o„,6„- (108) 

In order to extract the large time and distance asymptotics of (107) it is necessary to sum an 
infinite number of terms. In analogy to the finite temperature case we proceed as follows. 

• Consider all contributions at a given order in the functions K ab {q). 

• At each order we encounter two types of divergences: (i) infinite volume divergences that arise 
because in the thermodynamic limit the scattering states are normalized to delta functions. 
These singular contributions are ultimately compensated by corresponding divergences in the 
normalization (108) and are dealt with by an appropriate subtraction procedure, cf. [89] for the 
finite temperature case; (ii) "infrared" divergences, i.e. contributions that diverge as t — » oo 
or x — > oo. 

• We isolate the terms with the strongest infrared divergences at each order in K ab (q) and then 
sum up all these contributions, cf. [95] for a similar calculation in the finite-temperature case. 



4-.1. Finite Volume Form Factors for the TFIM 



As is clear from the above discussion the basic building blocks are the form factors. For the case of 
the transverse field Ising chain these naturally depend on the precise choice of basis of free fermion 
scattering states \k\, . . . , /c n )R,NS- For a particular such choice, the finite- volume form factors of the 
spin operators af have been determined in Rcfs [96, 97, 98]. The non- vanishing matrix elements 
are 

, 92>? \ P 1, • • • ,Pm) R = e-"K=i Pm] 

2 2n 

3 = 1 



m (qi, 



X I 



|n+f J 



(AJ 2 hY 



Le 



<i.i 



n 

i=i 



e Vp i 
Le 



pi 



>< n 



sin(^) 



-Hi, 



n 

Kl'=l 



sin(^ 



-ViPi> 



2n m 

nn 

j=iz=i 



sin 



where m is even (odd) for h < 1 (h > 1) and 

e k = 2J v / l + / i 2 -2/icosfc , e fc>fc , = £fc ± £fc ' , ^ = |1 - /i 2 ] 1 / 4 , 

n 



PER 



gGNS 

n e k, P 

p£R 



(109) 



(110) 



<j,<j'eNS p,p'eR 

For large L we have with exponential accuracy (in L) 

frfsl e^wl. (Ill) 

As we are interested in the thermodynamic limit we will set these terms equal to 1 in what follows 
(their deviations from 1 give rise to sublcading contributions). We stress that as a consequence of 
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the different quantization conditions of R and NS momenta there are no singularities in Eq. (109) as 
long as L is (large but) finite. The free fcrmionic basis used in Ref. [98] differs from the one discussed 
in Appendix A. However, the initial state can still be represented in terms of eigenstates of the post 
quench Hamiltonian H{h) using (B.8) and (B.10), provided that K(k) is chosen appropriately. 



4-2. Quench within the ferromagnetic phase: time evolution of the order parameter 



We now consider a quench within the ordered phase: we prepare the system in the ground state of 
the Hamiltonian H(ho) and at time t = suddenly change the magnetic field from ho to h, where 
ho, h < 1. We are interested in the case where the Z 2 symmetry of H{ho) is broken spontaneously 
in the ground state | \&o) - This is possible only in the thermodynamic limit. On the other hand, we 
would like to keep the length L of the system very large but finite in our calculations for technical 
reasons. In order to be able to work in a large, finite volume, we therefore take our initial state to 
be of the form 



|*o) = ^|{|0;Mns + |0;/io)r}, 



(112) 



where |0; /io)ns an d |0; ho) R are the ground states of H(h ) in the sectors with even/odd numbers 
of fermions respectively, see Appendix A. 4. For finite system size L, the state (112) is a particular 
linear combination of the ground state and first excited state of the Hamiltonian H(ho). On the 
other hand, in the thermodynamic limit L — > oo (112) becomes the symmetry broken ground state 
of Ho- As is shown in Appendix B the time-evolved initial state 

|*o(*)> = e- iff W*|* > (H3) 
can be expressed in terms of eigenstates of H(h) as 



where 



\B(t)) K 



l*o(*)> = 



e" 4E o* exp 



V2 



\B(t)) 



NS 



\B(t))j 



Vns{B\B) m y/ R (B\B) R 



[i ]T K{p)t 



-2is„tiA 



—pp 



\0;h) t 



R, NS. 



(114) 



(115) 



0<p£a 



For the choice of basis underlying (109) the function K(k) is given by 

*(*) = S[ f k) {h - ko) . (116) 

e ho (k)e h (k) (2 J) + 1 + hho - {h + h ) cos(fc) 

We note that this agrees with what would be obtained using the choice of fermions presented in 
(Appendix A). The expectation value of af is then given by 

n(B(t)\a*\B(t)) NS + NS (B(t)\a*\B(t)} R 



^Vr(B\B) r ns(B\B) ns 
We note that the diagonal contributions vanish 

a <£(t)k/W)>a = 0, a = R,NS, 

because 

e <jg e = — o g , 
where N is the fermion number operator (A. 6). 



(117) 



(118) 



(119) 
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The normalization is readily evaluated using the explicit representation (115) 
Ns( 5 l B >NS=ex P [ E log(l + ^ 2 (9))] , 

CKqeNS 

R (B|B) R = cxp[ ]T log(l + K 2 (p)) 
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CKpGR 



(120) 



so that, for large L, the norms rns(^I^)rns are approximately equal (we note that K(0) 
K(tt) = 0) 



/ f 71 ' dk \ 
r(B\B)k « ns(5|S) NS « cxp(i ^ — log(l + K 2 {k))) . 

= l + 0(e- aL ), 



(121) 



Importantly we have 

ns( b I b )ns 

where a is a positive constant. Eqn (122) allows us to simplify the expression (117) for the 1-point 
function to 

R (B(t)\a*\B(t)) NS 
r(B\B) r 

As we are interested in the thermodynamic limit this representation is most convenient and we will 
use it in the following. The normalization can be expanded in powers of K as 

1 2 

E ^ 2 ( fc ) 

_0<keR 



0(L- 2 ). 



(122) 



(123) 



,(B\B) R = 1+ J2 K 2 {k) + \ 



o<fceR 



-\ E K\k) + ... 



CKfceR 



i + E T2 « > 



(124) 



n=l 



where T 2n collects all the terms in which 2n functions K(k) are multiplied together. The following 
representation of T 2n turns out to be particularly useful 
i n 

r ^=^ e n^)- ( i25 ) 

'o<fci,...,fe„eR j=i 

Here indicates that the sum is only over terms with fc; ^ k m for Z 7^ m = 1, . . . , n. Eqn (124) is 
a formal series expansion as each term grows as L n , i.e. diverges in the thermodynamic limit. The 
divergences in the normalization (124) are mirrored by infinite volume divergences in the numerator 
of (117). Using (115) to expand Ng (.B(t)|o-£|.B(t)) R in powers of the function K{k) gives 



(126) 



Re[ NS (B(t)\tf\B(t)) n ] =2m x Al + R 1 {t)+ R 2 (t) + T 2 + R 3 (t) + T 2 R 1 (t) 



+ 



R 4 (t) + T 2 R 2 (t) + T 4 



where mg = (a^) /2 is the magnetization at the initial time and the R n (t) are finite in the 
thermodynamic limit. In terms of the "connected" contributions R n (t) the order parameter 
expectation value (123) is expressed as 



(*o(*)|of |*o(*)> = 2mg 



1 + E^i) 



n=l 



(127) 
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Eqn (127) constitutes a linked cluster expansion, where the contribution R n (t) is of order 0(K n ). 
As is shown in Appendix B (see Eqn (B.12)) physically the formal expansion in powers of K{q) 
corresponds to a low-density expansion, where the small parameter is the density of excitations of 
the post-quench Hamiltonian in the initial state. 

The next step is to determine the functions R n (t). Expanding the boundary states in 
NS <B(t)kf|S(t)) R we obtain 

oo ; n I 

m (B(t)\at\B(t)) R = ]T — Yl \{K{ qo )\{K( Pl )e^^-^^ 

l,n=0 ' ' fli.-,? n eNS j=l i=l 
Pl,...,pj€R 

x ns( _ 9i.9i,- •• ,-qn,q n We\Pi,-Pir ■■ ,Pi,-Pi)r 

oo 

= E C K0- ( 128 ) 

l,n=0 

In the following we consider the first few terms in the expansion (128). We will find that most terms 
exhibit long-time (infrared) divergences, the strongest of which occur in "diagonal" contributions 
n = I. We will determine and then sum these leading singularities to all orders in the expansion 
(128). 

4.2.1. Order O(K ) Contribution (n = I = 0). The zero-particle contribution in (128) is the 
ground state (of the post-quench Hamiltonian H(h)) expectation value 

C (0 |o) = 2mg = = (1 - h 2 ) 1 ^ . (129) 
4-2.2. Order O(K) Contributions. At first order in K there are two contributions 
C(H0)=i E K{q)e^ N s<-g,«k?|0> R = 4J%/^ E K{q) ^ q e^-, 

C (m = -i E K( P )e-^ NS (0\aZ\p,-p)K = U 2 h^ E K{p) j' mp . (130) 

0<pGR 0<p£R P 

In the limit L — > oo we can turn the momentum sums in (130) into integrals, which gives 

I{t) = ± hm [C (m +C m) ]=8J 2 h ff K{k) f nk cos(2ts k ). (131) 
For late times we can evaluate this integral by a stationary phase approximation 

»> + (132) 

where 

hJ 2 K'(k) , n 

Comparison of (128) and (126) then allows us to identify the function Ri(t) in the thermodynamic 
limit 

R 1 (t)=I(t). (134) 
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4-2.3. Order 0(K 2 ) Contributions. To order 0(K 2 ) there are two types of contributions. C( 2 |o) 
and C(o|2) are finite in the thermodynamic limit and well behaved at late times. They don't play a 
significant role in the following and we therefore refrain from presenting explicit expressions. The 
most important contribution to order 0{K 2 ) is given by 

C WD = E K(q)K(p)e- 2it ^-^ NS (-q,q\^\p,-p) R . (135) 



0<q€NS 
0<p£R 



The relevant form factor is, cf. eqn (109), 



NS\-9>9l°/b' -P>R - sin 9 sin p c ?p . 



where we have defined 

Cfcfe' 



I 



(136) 



(137) 



cos k — cos k' 

For large t and L the momentum sum in the NS-sector can be evaluated using Lemma 1 (D.l) and 
retaining only the leading terms in t and L. This results in 



C(i|i) = y/i E 
o<peR 



1 



4ie' 



+ ... 



(138) 



Crucially, this contribution exhibits an infinite volume divergence (for L — > oo) as well as an infrared 
divergence (for t — > oo). To 0(1) (in both L and i) (138) can be expressed as 



f w dk 

C ( i|i) = V£T 2 - V£ / -if 2 (fc)2te' fc + ..., 



(139) 



where T 2 is defined in (124). By comparing (139) to (126) we can determine the function R 2 (t) for 
large t 



f" <\k 

R 2 {t) = - / —K 2 (k)2te' k 
Jo n 



+ . . 



(140) 



where the dots denote terms that are subleading in t at late times. We note the leading term (140) 
arises from the region p rts q in (135). This observation will be important in what follows. 



4-2.4- Order C(if 3 ) Contributions. To this order there are several contributions. CV310) and C^o|3) 
are completely regular and do not play an important role in the following. The evaluation of the 
other contributions follows closely the calculation for C(i|i) and results in 

C\ m + h.c. = Re E iK(k 1 )K(k 2 )K(p)e 2it ^ k i +Ek 2-^ 

0<fei, 2 €NS CKpGR 

x Ns(-ki,k 1 ,-k 2 ,k 2 \o-^ l \p,-p) R 



R- c (C(o|i) 



T 2-^ e 4te ;^ 2 (p) 



0<p£R 



+ ... , 



C( 2 |i)+h.c. = Re E iK{k 1 )K{k 2 )K{p)e 2it ^ +£k ^ 

o<fei, 2 eR o<psns 

X R (-fci,fci,-fc 2 ,fc 2 |(7^|p, -p) N S 



= Re(C(i| )) 



T 2 -i 4te' p K 2 ( P ) 



+ ■ 



0<p£NS 



(141) 
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Comparison with (126) then gives 



Rait) = -I(t) f V d ^2te' v K\ v ) 
Jo w 



(142) 



where the dots indicate contributions that are sublcading in t. We note that (142) arises from the 
regions k\ t 2 ~pin (141). 



4-2.5. Order 0(K 4 ) Contribution C/212) - The most important contribution at order 0(K 4 ) is 
C(2|2) an d we now discuss its evaluation in some detail. From (128) we have 

K(q)K(q')K(p)K(p')e- 2U(£ " +£ r-'-^- £ ^ 



C, 



(2|2) 



2! 2 



o<<;^«'eNs 
o<p^p'eR 



NS\ 



q, q, -q', q'\o-'t\p, -p,p', ~p') R , 



(143) 



where the form factor is given by 



NS 



Each facte 



(-3, q, -q', q'Wi\p, -p,p', -p'h 



16\/£ ^q : p^q,p'^q' ,p^q' ,p' 
T4 F 2 F 2 f 2 F 2 4 ~A~ 

q q p p q,q p-.P 



x sing sing sin p sin p' 



2222 

C C ,C 1 C 1 1 
l^qp ^qp'^q'p^q'p' 

r 2 2 

qq' pp' 



(144) 



is associated with a singularity in the thermodynamic limit and it is useful to 



isolate these poles as functions of the NS-sector momenta, e.g. 

CpqCp'q — Cpp' {Cp'q ~ ^-pq) i 



2 2 2 
f r . — c , 
pq p q pp 

2 2 2 2 
^qp^qp'^q'p^q'p' 



Cp' q 2Cppi(c p i q C pq ) 



— C qp C q'p' C qp' C q'p ^ C pp' i C P' q c pq)( c p'i}' c pq' ) 



(145) 

-qq'^pp' 

Using the symmetry of (143) under exchange of p and p' the contribution of the last term in (145) 
can be reexpressed by substituting 



C qp C qp' C q'p C q'p> 
C 2 ,C 2 , 

qq'°pp' 



^ C qp C q'p' ^ C pp' C Pq C Pq' ^ C pp' C pq C p' q' ■ 



(146) 



The sums over the NS momenta can then be carried out using Lemmas 1 (D.l) and 5 (D.10). We 

[2 2] 

begin with the contribution due to the first term in (146), which we denote by C^y Here the 
superscript indicates the pole structure, namely two double poles. The NS sector momentum sums 
in (143) are performed by using (D.l) and retaining only the leading terms in L and t, which gives 



C 



[2,2] 
(2|2) 



0<p^p'£R 



4*4 n 



K\p') 



0<p£R 



1 - 



L 

8te' 



(147) 
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Working out the contributions arising from the other terms in (146) is more involved. Carrying out 
the sums over q and q' using (D.10) and retaining only the leading contributions in L and t gives 



4c; 



pp'CpqCpq 



0<p#p'£R £ P £ P' 

4 ^ K\p)K*(p')s\np^p'cl pl . 



Acp\p,CpqCp, ql K z (p)K z (p')smpsmp'c z pp , . (148) 

The leading contribution at late times and large L can then be extracted by using Lemma 2a (D.2) 
C 8i?) = -1T £ + • ( 14 9) 

0<p£R 

Here the superscript [1, 1] indicates that we are considering the contribution from two single poles. 
Putting everything together we arrive at the following result for Cr 2 \2) 

C(2|2) = egg' + eg = v^[k^(t)) 2 + r 2 i? 2 (t) + t 4 



(150) 



This allows us to identify the leading contribution to the term R^t) in (126) as 
R,(t)= l -(R 2 {t)f + .... 



(151) 



4-2.6. Exponentiation of the Contributions C(„|„). Our analysis of the first few orders in powers 
of K reveals the general structure of the expansion (128): at each order (except the very lowest 
ones) there are late time divergences that become stronger at higher orders. Moreover, the leading 
singularities are found in the "diagonal" contributions C(„i„) and we will now isolate these singular 
terms. At order 0(K 2n ) we have 



a 



(n\n) 



7.12 



ft 

]T' \\K( qj )K( P] )e^^-^ 



Ql,---,4n€NS j — 1 



X Ns( _ 9l. 91. ■ • ■ . -Qn,q n \o-t\pi,-pi ...,Pn, ~Pn) K , 

where the form factor in the limit of large L is given by 

Ns( _ 9l'9l' • • ->-Qn,Qn\^e\Pl,-Pl ■ ■ ■ ,Pn,-Pn) R = 



(152) 



4Vf 

L 2n 



n 



e 4 c 2 

Qj,Pl QjPl 



n sinfa)sin( Pj -) 
p2 P 2 



n 



i 



e 4 e 4 c 2 c 2 

j<j'=l Qj>Qj' Pj'P}' Ijlj' PiPji 



(153) 



The most singular terms at late times arise from the regions in momentum space 

Pj ss g P(i) , (154) 

where P is a permutation of 1, . . . ,n. There are n! such contributions with only double poles and 
they are all equal. They are determined by replacing 
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and then carrying out the NS-sector momentum sums using (D.l). This gives 

c [2,...,2] _ 1 4" V| 



(n\n) n \ £ri 



E n*o*) s 

0<p 1 ,...,p„eR i=l 



— - te' 
4 Pz 



+ ... . 
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(156) 



This is in agreement with (138) and (147). Using (125) we see that the ^-independent part of (156) 
equals T 2 „, so that 



E^) 2l | fc0 = ^ex P [E Iog(l + ^(p)) 



n=0 



+ ... . 



(157) 



0<pGR 

For t > we may invert the steps used to express (120) in terms of (125), which gives 

4tel 



£ C &N 21 = ^exp[£ log(l + *»(p) 



n=0 



1 - 



0<pGR 



+ ... 



V^ R <B|B> R expU f ^- 
L Jo ^ 1 



r7r dfc K 2 (k) 



1 --nm l2£ ' m 



+ ... . 



The contribution of all these terms to the 1-point function of the order parameter is thus 

r w dfc K 2 {k) 



7T l + if 2 (fc) 



|2 £ '(fc)| 



(158) 



(159) 



4-2.7. Exponentiation of the Contributions C\ n ±i\ n ). The leading terms (at late times and large 
L) in C( n ±i|„) can be summed to all orders in a similar way to our treatment of C( n \ n y The 
contributions C( n+ i| n ) are given by 



Re C( n+ i| n ) 



(n + l)!n 



0<9! ,...,«„, 9GNS 
0<P!,...,p„eR 



JjA-( % -)^(Pi) 



2 J te,+2iti:* , = 1 (e, I -£p s ) 



X Ns( _ 9l> 91, • • • , -«n, in, ~q, q\<?l\Pl, ~Pl, ■ ■ -Pn, ~Pn) R , 

where the form factors in the limit of large L are 

Ns(~gl» 91, ~Qn, Qn, ~q, gkf -pi, ...p n , ~Pn) R = -4iflJ 2 +1 JJ e* )Pj JJ £~^£~ 



(160) 



x n c 9Aj' e p*pj' n sin ^) sin ^)fr 



nn 2 

j,i=l c g 3 pj TT ^g, Pj am (j -pr 

n" r 2 „2 11 ,4 F 2 11 „2 

i<j' = l ^Ijli'^VjVji ,-=! C 9>9j C <? 7 = 1 1 



£2r; 

j',Z=l m=l 
n .4 . n 2 

£ <?,p 3 sm 9 tt c m 

,2 



(161) 



i=l "J^J —J- >«^y j- 

Following the same reasoning that led to (155), we focus on the contributions arising from the 
regions 

Pj ~ Qj i i = !,••• , (162) 

where (qi, . . . , <? n +i) is an arbitrary permutation of (q\, . . . , q n , q). All of these contributions are 
the same and hence the leading behaviour at late times can be extracted by the substitution 

n-=inr=i^ 



[" TT™ r 2 n r 2 

lj=l 111=1 ^qjVi TT IPi 



^■j<j' = l C %<lj' C PiP 3 ' 7 = 1 C; 



n^^(-+i)'ii c L i +--- 



(163) 



HI 



i=l 
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The NS-sector momentum sums can be carried out using (D.l), which gives 
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K<3 C (n+l\n) 



1 4™+ 1 W 2 v / C 



L n+1 



^)cos(2te g )^n^) 



o<p 1 ,...,p n eR 



9 i=l 



4 - 



= -4W 2 ^ 21 ^gif(fc)cos(2 tefc ) Sin/ ' 



Summing over n then results in 

Re E C fcifi) - V^r^Wr exp 



n=0 



dfc if 2 (fc) 

7l + if 2 (J;) 



\2e'(k)\ 



+ 



(164) 



(165) 



where I(t) is given by (132). The analysis for the contributions C(„|„ +1 ) is completely analogous 
and finally leads to the following result 



Re£ C( " +1 ';i + D C( "'" +1) = yTW) exp 



n=0 



<S|B> R 



T dfc K 2 {k) 
+ K 2 {k) 



r dfc 



|2e'(t)| +..(166) 



We see that by virtue of the extra factor I(t) this is always subleading compared to the "diagonal" 
contribution (159). 

4-2.8. Form Factor Result for the 1-point Function (^o(t)Wi\^o(t)) Combining the results in the 
previous subsections we conclude that 

(*o(*)|of |*o(*)> ~ [1 + I(t) + Oit- 1 )] e"*/r , (167) 

where the decay time r is given by 

r* dfc 



r = 



-K 2 (k) [2e'(k)} + 0(K 6 ). 



(168) 



Here we have used (151) to infer the absence of 0(K 4 ) corrections to the decay time. The corrections 
in the prefactor are expected to be as there are such (subleading) contributions in the various 

terms C( n |j) we have calculated above. We note that these terms are larger than I(t) at late times 
as I(t) ~ i- 3 / 2 . 

The result (168) is in agreement with the one obtained by applying the cluster decomposition 
principle on the asymptotic result for the two-point function obtained in the determinant approach 
(cf. Eq. (69)), which gives 

l + if 2 (fc) 
1 -K 2 {k) 

Indeed, expanding (169) in the limit where K 2 (k) -C 1 gives precisely (168). We note that (168) 
is a rather good approximation for a wide range of magnetic fields h a and h. For example, for a 
quench from ho = 0.1 to h — 0.8 the relative error is less than 1%. This shows the effectiveness of 
the form-factor approach. 



J0 7T 



s'(k). 



(169) 
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4-3. Two-Point Function in the Ordered Phase 

We now turn to the time evolution of the two-point function for a quench within the ordered phase. 
The quantity we want to evaluate is 

p**(t,l) = (Mt)\< +( or x m \Mt)) 



+ 



;ns 



(170) 



(171) 



2 R {B\B) R 2 NS (B\B) NS 

In the large- L limit both terms contribute equally, so that 

lim ^( M )=lim »<*('>l^ n f !*('»». 

L^oo L^oo R (B\B) R 

We note that as a consequence of reflection symmetry the sign of £ does not matter, i.e. 
p xx (t, £) — p xx (t, —£). We therefore will assume from now on that 

£>0. (172) 

The Lehmann representation for the numerator on the right hand side of eqn (171) is 



R (B(t)\a x m+t a x m \B(t)) R = £ 



v i n ~ l 



1 



n k ^ 



3=1 



Y[K( Pi 



n\ l\ ' — ' * — ' r! 

I n =0 o<fe 1 ,...,fe n 6R o lr .., 9r £NS 

o<pi,...,pj eR 

x e -2ftEj=i ^-E: =1 £ M R (-fc 1; fc l7 . . . , ~k ni k n \a x m+l \qi, • ■ ■ , 9r>NS 
X NS (g r , . . . ,<7i|cr^|pi,-pi . . . ,pt,-pi)n 

OO 

= X] D (2n\r\2l)- (173) 
Z,n,r— 

Like in the case of the 1-point function we focus on the terms with the strongest possible singularities 
in the form factors for a given order in the formal expansion in powers of K(p). These are obtained 
by taking 2n = 21 = r. 

4-3.1. Contributions at O(K ). These are equal to the ground state two-point function 

r(0K + ,<|0)r, . (174) 
Using a Lehmann representation to evaluate (174) gives rise to contributions 

1 



D 



(0|2s|0) 



E 



0<5i,...,g 2a eNS 



(2«)1 



R<°l CT m+<ki. ■ • ■ ,<?2 S >NS ns(<72 S , ■ • . ,(?l|<JO) R . (175) 



For s = we have 

^(OIOIO) = 
while for s = 1 we obtain 



(2mg) 2 =e=(l-^ 2 ) 1/4 , 



(176) 



97,2 jit fir I 9i-g2 \ 

D m = / d Ql dq 2 [ I > e«(*+*> + O(L^) oc h*. (177) 

^ ^-TT £ qi £ q2 e q 1 ,q 2 

We see that at large £ £>(o|2|o) is exponentially small compared to -D(o|o|o)- The contributions 
-D(0|2s|0) with s > 1 are suppressed by additional exponential factors. We conclude that for large 
distances £ we only need to retain -D(o|o| o) - 
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4-3.2. Contributions at 0(K 2 ). At order 0(K 2 ) the largest contribution at large distances and 
late times arises from 



D 



(2|2|2) 



= E o E ^(fc)^(p)e- 2lt[£ ^ £fcI e" (9l+92) R,(-MKki,g 2 )NS 



0<fc,pGR <7i,<j 2 eNS 



x Ns(?2,gikmb.-p>R . = d(1) 



(2|2|2) 



(2) 

(2|2|2)' 



(178) 



where £> 



(1,2) 
(2|2|2) 



denote the contributions arising from terms with k ^ p and k = p respectively. Using 



(D.4) to carry out the summations over q\p we arrive at 



D 



(i) _^ 

(2|2|2) L 2 



E 



-2it [Sp -e k] (Sk+£ P ) 2 K(k)K(p) 



CKfc^pGR 

+£k£ P 



(cos k — cosp) 2 
'k 



sin k sin p e\ 



sin 2 cos ^(jfe + p) ) _ sin 2 cog _ p) ) 



+ ■ 



(179) 



We note that in this expression there is no singularity if we consider the limit k—p — > 0. Nevertheless 
for large I and t the leading contribution to the double sum arises from the vicinity of k = p. In 
order to isolate this contribution we turn the sum over k into an integral using the Euler-Maclaurin 
sum formula and then deform the integration contour into the upper half plane. As e' p > for all 
p > the contribution of the first term in square brackets will be negligible for large t. The same 
holds true for the second contribution as long as 2e' p t ± I > 0. On the other hand, if 2e' p t ± £ < 
we need to deform the contour into the lower half-plane (in the variable k). In doing so acquire a 
contribution from the double pole at p = k. The residue is dominated by the factors involving t 
and I as both of these are assumed to be large, so that we end up with 

4£ 
L 



^(2|2|2) 



E K 2 (p)e H (e-2e'( P )t)[2e' p t -£]+... , 

0<pGR 



(180) 



where H {x) denotes the Heaviside step function. The second contribution £v 2 |2|2) * s gi yen by 



D 



(2) 

(2|2|2) L 4 



8£ sin 2 (A:) K 2 (k) 

^ ^Hk) ^ 



(91+92) sin 



in 2 (^) 



4 4 
e q 1 ,k € q 2 ,k 



2 2 
C kq 1 C kq 2 - 



(181) 



CKfceR ' 9i, 2 €NS qi q2 «!^ 2 

Using (D.6) to carry out the summations over qi^ we arrive at 

l(2) 



(2|2|2) 



£ E !"£ +••• ( 182 ) 

where we again have only retained the leading terms at large L and I. Adding the two contributions 
(180) and (182) we obtain 



£>(2|2|2)(M)=£ 



T 2-T E ^ 2 (fc)(^(^-2te' fe )[2te' fc -£]+. 



0<fceR 



+ ... 



T 2 - y E R2 ( k ) ( 2t£ k - 2te'{k)) + 1 8 H (2te' k - £)) 



fc>0 



(183) 



This contains a part T 2 that diverges in the thermodynamic limit, but which will again be 
compensated by the denominator in (171). 



Quantum Quench in the Transverse Field Ising chain I 
44. Exponentiation of D {2n \2n\2n) 

The dominant contribution at order 0(K 2n ) is again the "diagonal" one 
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D(2n\2n\2n) — 



z' n 



(n!) 2 (2n)! L 4 « ■ 

V / V / 0<fci,...,fc„eR j-=l 
0<pi Pn€R 



sin(fc r ).ft'(A; r )e ^ sin(p r )i ; ir(p r )e" 



-2ite v 



n 



n 



1 



4c? 



4c? 



;</ , fcjfc,, fc,,fc,, 'PiPj/ Pi,p,/ g 1) ...,q 2 n£NS 
2 ™ sm 2 / gl^gj/ \ " 2n 

x n J 2 nn 4e i.«. c Pr,g. c fcrg.<w 

«i' 9;' r =ls=l 

(184) 

In order to extract the large time and distance behaviour of (184) we need to focus on the regions 
where the form factors exhibit the strongest singularities. For the case n = 1 discussed in subsection 
4.3.2 these regions were 

• qi w ±k , <72 ~ TP and P ~ k , 

• g 2 ~ ±fc , (71 w =FP and p w fc. (185) 
For general n we need to focus on the regions 



Pi 



/ qsi 



\qs 



-o-\Pi 



0~nkQ n 

V -o n p n 1 



(186) 



where Uj — ± and (Qi, . . . , Q„), (Si, . . . , Sin) are permutations of (1,2 ... ,n) and (1,2..., 2n) 
respectively. By symmetry all these regions contribute equally, which gives rise to a combinatorial 
factor of n!(2n)!. We therefore focus on the single region where (Qi, ■ ■ ■ ,Q n ) = (l,---, n ) an d 
(Si, . . . , S2n) = (!)•••; 2n), together will all regions obtained by exchanging any pair S2 r -i <^ 5*2r- 
We start by expressing (184) in the form 



D 



(2n\2n\2n) 



L 4n n\2 r ' 



E' e 2it ^n=i(=*,-=p,) f(ki,...,k n ;p 1 ,...p n \q 1 ,...,q 2 n) 



0<k 1 ,...,fc n £R 
0<P!,...,p n ER 



gi,...,92neNS 



n p i(q2 r -i+q2r)e 
e t -fer92r-l t -fer92r C Pr92r-l L -Prg2r' 



(187) 



where it is understood that the summation only terms arising from the region we consider are 
retained. We then carry out the sums over the q s 's using (D.5) and (D.9). When doing this we 
need to distinguish the cases k r ^ p r and k r — p r as they give rise to different kinds of singularities in 
the q s summations. For k r = p r the summation over both qi r -i and qi T gives a leading contribution 



C/2r-l,g2rGNS 



p r —k r 



L 
16 



1 - 



f(Pr, ~Pr) + f(~Pr,Pr) 



sin (p r ) 



+ ... 



(188) 



k r —p r 
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where we have defined 

f(P,0) = f(h, ■ ■ -,k„;pi, . . . ,Pnkl, ■ ■ .,q 2 r-2,P,q,q2r+l ■ ■ • ,92n)- (189) 

The extra factor 1/L in (188) is present as we are considering a single term in the sum over k r . On 
the other hand, for p r ^ k r we obtain 



}_ e 2it[e kr -e Pr 
r ,3 Z—d 



I? 



0<fc r SR 



E 

<?2r-l,92r6NS 



i(q 2r - i+q2r)i ri r , r r 



: E 

o<t r eR 



2it[e kr -e Pr ] 2 



p r fc r 



x /(fci, . . .,fe„;pi,.. . • • • ,?2n) 

4 sin £; r sinp r 



4 sin (p r ) 



(190) 



Similarly to our treatment of the n = 1 case in subsection 4.3.2 we can now carry out the summation 
over k r with the result 



10 0+„' UO += .' mf(PT,-Vr) + f(-Pr,Pr) 

e H (i-2ts Pr )[2ts Pr -i] 4sh?(p ^ 



+ ... . 



(191) 



k r —p r 



Combining the two contributions (188) and (191) we find that the leading terms at large I and t 
are 



f(Pr,-Pr) + f(-Pr,Pr) 



— 'fir 



Carrying out all q s and ki summations in the same way gives 

£ n 1 

I>(2„|2„|2n) ^ E II [ L - 4£ - MhV 2ts' pT )[2te' Pr - £}] T ^— i 



+ ... . (192) 



(Pr) 



0<pi,...,p n eRr=l 

X E /(Pl>'--P™;Pl J ---Pnl Cr lPl:- Cr lPl ; --- J Cr ™Pn ; - (7 nPr l ) + --- 



cri,...,(Tl=± 



n! 



E' U k2 m 



it 4 

i -- - ^e H (e-2ts' Pr )[2te' Pr -e} 



(193) 



0<pi,...,p„eRr=l 

Like for the 1-point function the sum over n can be taken by inverting the steps used to express 
(120) in terms of (125), which gives 



OO / 

E D (2n\2n\2n) = £ exp I ^ In 



n=0 



v 0<peR 

- £ ns( b I b )ns ex P 



1 + X 2 (p) ( 1 - - (I + 9 H (£ - 2te' p )[2te l p - £]) 



+ ... 



■ J2 K 2 (p)(e + e H (£-2ts' p )[2ts' p -£]) 

0<pGR 



(194) 

Here we have retained only the G(K 2 ) term in the exponent as the higher orders are beyond the 
accuracy of our calculation. This then gives the desired result for the two-point function (171) in 
the ordered phase 



lim p xx (t,£) ~ (l-h 2 )* exp 

L— >oo 



2 r^K 2 (p) {£ + e H (£-2te' p )[2te' p -£]) 
Jo 71 



(195) 
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Eqn (195), including the prefactor, is expected to be accurate at late times and large distances 
as long as we do not quench too close to the critical point. 

4-5. Two-Point Function for Quenches in the Disordered Phase h ,h > 1 

As a consequence of the Z 2 symmetry the 1-point function of crj is identically zero for quenches 
within the disordered phase, i.e. ho, h > 1. We therefore turn to the two-point function. As shown 
in Appendix A the ground state in the paramagnetic phase is the NS vacuum |0; /io)ns- Hence the 
2-point function after the quench is equal to 

pII(M) = N S (%|» i (196) 

In the basis underlying the expression (109) for the form factors, the boundary states are again 
given by (115), but now the function K(k) is given by 

K(k) = Sl ; (fc) {H ho) . (197) 

e ho (k)e h (k)(2j) +l + hh -(h + h ) cos(fc) 

The Lehmann representation of the numerator is again given (173) if we replace R — > NS. The 
difference to the ordered phase is that now only form factors with odd numbers r of particles in the 
intermediate states are non-vanishing. 

4-5.1. Order 0(K ) Contributions. The 0(K ) contribution obtained upon expanding the 
boundary states in powers of K is 

NS (0K + ,(i)<(t)|0) NS - ns(0K +£ <J0) ns . (198) 

This is equal to the static zero temperature correlator in equilibrium. Inserting a resolution of the 
identity we conclude that the large-^ behaviour in the L — > 00 limit is determined by one-particle 
intermediate states 

1 e* 9 '^' C dq e iqi 



£(0|i|0)(M) = t^Ph- ]T — - = 23 A \ O(L-i) 

(199) 

where we have defined 

A = tVh. (200) 

The integral can be carried out approximately by bending the contour around the branch cut of 
the energy in the upper half plane, which gives 




D mo)«^I^Tj h -JT = v{i) - (201) 

We conclude that the leading contribution in K{q) is exponentially small, which means that we need 
to evaluate all sublcading terms with the same exponential accuracy as well. The contributions due 
to 3,5,. . . particle intermediate states are small for large I when compared to (201) 

^(0| 2s+ i|0)«^ (2s+1)W , (202) 
and can hence be ignored for our purposes. 
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4-5.2. O(K) Contributions. The leading contributions to order K are obtained by taking single- 
particle intermediate states into account only. This gives 



D. 



(2|1|0) 



+ D 



(0|1|2) 



— J2 -r 1 sm(2te k -k£) + .... 



(203) 



The corrections to the rhs of (203) can be evaluated using contour techniques but are negligible at 
large distances. As a function of t for fixed n -D(2|i|o) + -C(o|i|2) displays oscillatory behaviour on 
top of a slow t -3 / 2 power-law decay in time. This can be seen by turning the sum in (203) into an 
integral and evaluating the latter by means of a saddle-point approximation, which gives 



D 



. j~. 2JA^ K(k a ) J 

(2|1|0) + £(01112) * —FT £ Im 



o=± 



£ k a 



-ik r J: 



where 



k± = arccos 



c 2 (l + h 2 ) + h 2 i 



h 



x = 



AJt 



In the limit I ft — > this can be simplified further 



-°(2|1|0) + D(0\l\2) 



£ 



(2Jt) 

£ 



3/2 



h- 1 



1/4 



h-h 



Ah(h - 1)V5F V 



sin 4Jt(h - 1) 



if 



(2Jt) 



3/2 



h-l 
h + l 



1/4 



h — h 



4h(h a + l)V5F V 



sin AJt(h + 1) 



AhJ t 

h+l £ 2 
AhJ J 



(204) 



(205) 



(206) 



4-5.3. 0{K 2 ) contributions. The next contributions we need to consider are second order in K(k) 
and arise from 

£ K(k)K(p) NS {-k,k\a^ +e (t)<(t)\p,-pW (207) 

0<fc,p6NS 

Inserting a resolution of the identity between the spin operators we see that for large £ and t the 
dominant contributions are generated by 1-particle and 3-particle intermediate states. The former 
is given by 

£(21112)= £ £*W*(p) NS<-MK-n(t)|g> R 

k,p>0 q 

x R (q\a^{t)\p,-p) NS . (208) 

Like in our analysis of the 2-point function for quenches within the ordered phase we again have to 
consider the two cases k ^ p and k — p separately. Denoting the corresponding contributions by 
respectively, we find in the limit of large £, t and L — > oo 



and n (2) 

^(2|l|2) dllU ^(2|1|2) 



D 



(2) 

2|1|2) 



2JA \ — ^ T , 9/ ,,e ifc l^l 



^(2|l|2) 



L 

8JA 

IF 

+ ... 



£ ^ 2 ( fc ) 



feeNS 



(209) 



E 

0<fc/p6NS 



K{k)K{p) 2 sin(fc) sin(p) 
cos(fc) — cos(p) 



2 C 2 



eie 



£fe sin(fc^) e p sm{p£) 
sin(fc) sin(p) 



2it[s k -e p ] 



(210) 
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In the thermodynamic limit (210) can be written as a double integral, but we have not succeeded 



in simplifying it in a useful way. The contribution D 



(2) 
(2|1|2) 



is time independent and exponentially 



small in £. In contrast to this, -D( 2 |i|2) displays power-law decay in t for fixed I 



D 



(i) 



oc t 



'(2|1|2) ^ L ■ ( 211 ) 

We now turn to the 0(K 2 ) contribution involving 3-particle intermediate states 

£>(2|3|2) = \ E E K{k)K(p) NS (-k,k\<7* n+e {t)\ qi ,q 2 ,q 3 ) R 



0<fc,peNS 9i,92,93£R 

x R(q3,q2,qiW^ n (t)\p,-p) NS . 



We denote the contributions from k ^ p and k = p by -D^i^) 
lengthy calculations we find 

4£ 



and D 



(2) 

(2|3|2) 



(212) 

respectively. After some 



D 



(2) 
(2|3|2) 



1 



o<fceNS 

4JA K 2 (k) „ i2 ^ 1 ^ AJ 2 K 2 (k) 

+ ~r E — ^cos(^)-4^ 2 -i)v(i)- E — j- Li + --- 



L 0<feGNS £k ^ 0<k€NS ° fc 

This is time independent and contains a piece that diverges with the volume as expected. 



(213) 



The contribution D 



(i) 

(2|3|2) 



is given by 



^(2|3|2) 



2JA 



COS (2t[£fe - £ p ]) 



0<fe#pGNS 



6 



E 



91,92, 93 eR 



n 

3<l 



2 /9j-9l ' 



9j,9l 



3 2 2 p iq m n 

J J ffe,9rr. ' /'->/. 



4c, 



m— 1 



(214) 



Carrying out the sums over g., using the same techniques as for the other contributions we eventually 
arrive at 



D 



(i) 

(2|3|2) 



8JA 



K{k)K{p) 2 sin(fc) sin(f>) 



£ 2 £ 2 



c fe,p 



cos(k) — cos(p) 



£fc sin(fc^) e p sm(pt) 



sin(fc) 



sin(p) 



„2ii[e fc -£„ 



0<fc^pGNS 

+... (215) 

The large time and distance behaviour is thus the same as for -D(2| i|2) - 

In summary, the combined 0(K 2 ) contributions can be divided into two categories: 

(i) Time-independent contributions 

dk K 2 (k) 



AJA 



I 

J —IT 



2tt £fc 



+ D. 



(0|1|0) 



■At 



^ dk 4J 2 K 2 (k) 
/ 27 e 2 (k) ■ 
(ii) Time-dependent oscillatory contributions 



4JA 



4(/i 2 - l)V{i) 

epende: 
r dk_ 
Jo * 



dk dp K(k)K{p) 2 sin(/e) sin(p) 
2 e 2 (k)e 2 (p) fc,p cos(fc) — cos(p) 



£fe sin(fc^) s p sm(p£) 



sin(fc) 



sin(p) 



J2it[e k ~e p ] 



(216) 



(217) 
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4-5.4- 0(K 3 ) contributions. In order to infer the structure of higher-order oscillatory contributions 
we determine the following 0(K 3 ) term 

JA^ T W \ ,-,/, N -rrl N Sffi k\ SU1 k% SU1 J) 3 2 

~pr~J2~J2 9k 1 ,k 2 9k 1 ,-k 2 



D iMm + D (2 | 3 |4) = ^ KikJK^Kip) 



91,2,3E R 



o 2 o 2 o 2 3 

^91,92^91,93^92 ,93 2it[e fel +£ fc2 -£p] TT q,2 



£tfi ^Oo £fl 



e^HK^v > ( 218 ) 



£=1 



where we have dehncd 



2 sin 



(^) 



£p + £fc 



(219) 



(%) Contributions for p 7^ fci,2- We first consider only the contributions with p ^ fci^ and denote 
them by -D^j^) + -^(2?3|4) ■ ^ n or d er to carry out the sums over qj it is useful to rewrite the product 
of pole factors XIi=i 11.7=1 ^ c h<u usm g the identity 

Ckiq c k 2 qGpq c k 2 ki ^ c pki [ c kiq c pq\ 

The fully decomposed expression reads 



Cpk 2 [ c k 2 q c pq] }• 



(220) 



n%Ck iq] c k2q] c pq] = [8c fe2fel ] 3 I -c 3 fci 



2^&i)fci>f ) ^fci,/ci,fci t^pp.p.ki H~ Cp,p,p 



where we have defined 



(221) 



(222) 



We then can carry out the sums over gj using Lemma 4 (D.6) and retaining only the pole 
contributions. 

A. The combined contributions to (218) arising from c 3 ki C kukuP , c 3 k2 C k2M:P , c^ ki c pk2 C kuP , kl and 
c pk 2 c pki^k 2 .p,k 2 hi the decomposition (221) are of the form 

D^IL + = ^ £ K( kl )K{k 2 )K(p) e ^(k l)+ e(k 2) -s p] 



'(4|3|2) ^(2|3|4) - L 3 



o<fei, 2 ^peNS 



e fei,p e fc 2 , P c fc2fci sm fci sinfc 2 



e fci e fc 2 £ p 



c pfel - c pfc2 ] sin + h.c. + . . . . (223) 



B. The contributions due to <? pk C-p,pM and <? pk £p,p,k 2 in (221) arc 



^(1,2) ,£,(1,2) 



8JA 



'(4|3|2) ~(2|3|4) - L 3 



E 



Kjk^Kj^Kjp) 2 «t[ e (fc 1 )+ e (fc 2 )- s „l e fci,P £ fc2,P 



Cfc 2 fei 



0<fci, 2 ^p€NS fel fc2 p fcl >' C2 

x sinp [sin(£/ci) sinfc2£fc 1 efe 2 , p Cp fel - {fci O fc 2 }] + h.c. + . . . . 



(224) 
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^(1,3) ,£,(1,3) 



IQJA 



'(41312)^^(21314)- L 3 



E 



K{k 1 )K{k 2 )K{p) 2»t[e(fc 1 )+e(fc 2 )-e p ] ffci i pffc2 i P 



0<k 



sin(^fci) sin/c 2 -^c„fc 2 - {fci o fc 2 } 



c fei,fe 2 
+ h.c. + ... 



x c fc2fcl sinp 

D. The contributions of the c 2 fci c p fc 2 C PjPj fc 2 and c pkl c pk2 C PtPtkl terms are 

K{k 1 )K{k 2 )K{p) 2 «t[ e (fc 1 )+ e (fc 2 )- e „l £ Lp £ Lp 



(225) 



n (i,4) n (i,4) _ 8J.A 

^(4|3|2) + ^(2|3|4) - L 3 ^ £ £ £ 2 

0<fe 12 # P eNS E ^k^ p 



E 



c fcl,fe 2 



x c fe2fel snip 



e 2 c 2 

sin^feiJsinfe-^-^-^i^M 



£fc 2 Cpfcj 

E. The contributions of the c 2 pki c pk2 C kl MM an d c pkl c 2 pk2 C k2M . kl terms arc 

\^ K{k 1 )K{k 2 )K(p) 2<t[ e (fci)+ e (fc 3 )- ep ] e Lg e Lp 

=-2 e ^2 



+ h.c. + ... . (226) 



n (i, 5 ) n (i,s) _ 8J71 

^(4|3|2) + ^(2|3|4) - L 3 



0<fc 



C p fe 1 C p fc 2 . 

x — — - — smp 

Cfeifc 2 



li2 ^ peN s £ ^ £ ^ £ v 



sin(£fc2)sinfci-^c p fe 1 - {fci o fc 2 } 



: fei,fe 2 

+ h.c. + . . . (227) 



F. Finally, the contributions of the <? pkl c pk2 C kl . k2 . p and c pk2 c pkl C k2 , kl , p terms arc 



^(1,6) ,£(1,6) _ 



32 J A 



E 



K(k 1 )K(k 2 )K{p) 



cos (2t[e kl +e k2 - e p ]j 



'(4|3|2) --(2|3|4) - £3 £ fci e fe2 £ p 

fei,2^p£NS fcl fe2 p 

x e fe 1 ,p e fc 2 ,p c fcip c fe 2 pSfei,fc25fc 1 ,pfffc 2 ,p sm (^(ki + k 2 +p)), 
where here the momentum sums are over the entire Brillouin zone. 



(228) 



The contributions from all other terms are subleading. We observe that the leading terms at large 
£ and t in (224), (225), (226) and (227) combine to 

Kik^Kik^Kip) 

'(4|3|2) T^ (2 |3| 4) ~ - - 



E 



a=2 



cos 



(2t[s kl 



< >P 4,P sinA: 2 ship 



c P k 2 [c pk2 - CpfcJ sin(^ci) + . . . . 



(229) 



fei,fc 2 



Both (229) and (228) can be simplified further, because for large t and i the dominant contributions 
arise from the "double pole" factors and can be extracted using Lemmas 2a and 2b of Appendix D. 
This leaves us with 

+ - ~ E E ™ - - ^ 2 (p) 



x 2 



0<pGNS feGNS 

[f+(2te;-^(2fe;-f)] 



+ ..., 



(230) 
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(231) 



Combining (230) and (231) we arrive at our final result for the leading asymptotics of the p ^ k\ i2 
contributions to (218) 



D 



(i) 



+ D 



(i) 



(4|3|2) r ^(2|3|4) 



AJA K(k) . , 

2^ — ^ sin (2te k - 



L z — ' £fr 
+ ... 

[■D(2|l|0) +-D(0|1|2)] 



^ K 2 (p)[2te' p -e]9 H (£-2te / p ) 



CKpGNS 



■j E * 2 (p) [2te' p -e]0 H (e-2te' p ) 



0<pGNS 



+ ... . (232) 



(ii) Contributions with p = k\ or p = k 2 . Using the symmetry under k\ -H- k 2 we can express these 
contributions in the form 



D {2) +D {2) 

^(2|3|4) + ^(4|3|2) 



0<fei, 2 GNS gi,2,3£R 

x NS (-k 1 ,k 1 ,-k 2 ,k 2 \a? n+( ,\q 1 ,q2,q3)R R.(<73, 92, qi KJ&1, - fc i>NS + h.c. (233) 
In order to carry out the qj sums we rewrite the pole factors as 

3 

II 8c *ntt c Mj = (8c fc2fcl ) 3 {c 2 (A:i,fci,fci) - c k2kl b{ki,k 2 \qi)cl iq2 cl iq3 - c^cl^bih, k 2 \q 2 )c 2 kiq 
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- Ck 2 k 1 cl iqi cl uq2 b(ki,k 2 , \q 3 ) + c 2 klk2 cl iqi b{k u k 2 \q 2 )b{k u k 2 \q3) 

+ c l 1 k 2 c t iq2 b ( k i,k 2 \q 1 )b(k 1 ,k 2 \q 3 ) + 4 ife2 4 ig3 6(A;i, fc 2 |gi)6(fci, fc 2 |g 2 ) 

-^ 2fcl K*i,*2|?i)6(fci,fe|«2)6(«:i,fe2|«3)}, (234) 

where we have defined 

b(kx,k 2 \q) = c kiq - c k2q . (235) 

The various terms in (234) contribute in qualitatively different ways, depending on their structure 
when viewed as functions of (the complex variables) qi, q 2 and 93. 

1. Terms with only simple poles give rise to contributions of order C(i _1 ) and can be ignored. 

2. Terms with only double poles arise from the first contribution on the r.h.s. of (234). Using 
(D.6) to carry out the sums over qj we obtain after some calculations 



D 



(2,1) 
(4|3|2) 



D 



(2,1) 
(2|3|4) 



1 

L 



V(£) J2 ^ 2 (fci)^(fc2)cos(2te fe2 )c fe2 , fcl 

0<ki, 2 €NS 



L i — ' e k 

keNS K 

which can be simplified further if required. 



4JA K 3 (k) . , n 

2^ / ' sin(2te fc -tk) + ..., 



(236) 
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3. Terms with only one double pole arise from the fifth, sixth and seventh contributions on the 
r.h.s. of (234). Using (D.6) and (D.4) to carry out the sums over qj we obtain after some 
calculations 

lUA J2 K 2 (k 1 )K(k 2 ) cos (2te k2 ) sin(fc 2 ) cos(M) 



D 



(2,2) 
(4|3|2) 



D 



(2,2) 
(2|3|4) 



0<fci^fe 2 GNS 



x c k2kl 



Here we can carry out the k\ sum using Lemma 3 D.4, which gives 



D 



(2,2) 
(4|3|2) 



+ D 



(2,2) 
(2|3|4) 



4JA K 3 (k) . . 
— -j^ 1 sin (2te fc -£k)+... . 



(237) 



(238) 



feeNS 



4. Terms with two double poles arise from the second, third and fourth contributions on the r.h.s. 
of (234). Using (D.6) and (D.4) to carry out the sums over qj we obtain after some calculations 



D 



(2,3) 
(4|3|2) 



D 



(2,3) 
(2|3|4) 



4JA 



E 



L kt^s £k 



K{k2) sin(2te(fc 2 ) 



E K2 ^) 

0<fci 6NS 



At 



(239) 



This is in fact the leading contribution. 

Combining (232) and (239) we arrive at the following result for our 0(K 3 ) contributions 



D 



(4|3|2) 



+ D 



(2|3|4) 



u)+£>(0|i|2)] \ E fj- £ - [2te'(p)-^]0if(^-2te'(p)) 



0<p6NS 



8J.4 K 3 (k) . , ni , 
+ — 22 — sin ( 2te fe -tk)+... 



(240) 



We see that as expected there is a contribution that diverges in the thermodynamic limit. However, 
in addition there are terms that become very large for I » 1 



4.5.5. Resummation of Leading £>( 2n +2|2n+i|2n) + ^(2n|2n+i|2n+2) Contributions. 

The above calculation of the 0(K 3 ) contribution shows that the dominant contributions arise from 

the regions 

/ o~ik Ql 



kQ! ~P 




-0-1P 

\°~2kQ 2 



O-j = ±, 



(241) 



where (Qi,Q2) and (Si, S2, S3) are permutations of (1,2) and (1,2,3) respectively. Motivated 
by this observation we therefore consider the analogous regions for the 0(K 2n+1 ) contributions 
-D(2n+2|2n+i|2n) and -D( 2n |2„+i|2 n +2) • They are 





/ 951 \ 

Qs 2 



QS 2n -i 

Qs 2n 



( °i k Qi \ 
-V1P1 

-O-nPn 

\a n+ ikQ n+1 1 



(242) 
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where <jj = ± and (Qi, . . . , Q n+ i), (Si, . . . , S^n+i) are permutations of (l,2...,n + l) and 
(1, 2 . . . , 2n + 1) respectively. Our goal is to determine the contribution of the regions (242) to 



D 



(2n+2|2n+l|2n) 



+ h.c. = C Re 



fcl ,...,fc n _ ( _ 1 GNS 
Pl,...,p„£NS 



n+1 

n 

.r=l 



sin(/c r )ii'(fc r )e 



2ii£fc r 



n 

.s=l 



sin(p s ).ftr(p s )e 



-2 '(/,£,, 



n+1 



n 



x n c 2 e 4 n c 2 e 4 e 

;</' ^hh< klM' m<rn' ^P™.P m ^P m ,P m < qi ,...,g 2n £NS L r=l 
2n+l 



'2n+l 



2n+l 



x n 2e L +1 , 9s c ^ + « 



]^[ ^ e k r ,q s e p r ,q 3 C k r q s Cp 



.r=l 



where the constant C is 
C 



AJA 



(243) 



(244) 



L 4n + 2 (n+ l)!n!(2n + 1)!' 

As all regions (242) contribute equally, we focus on the case (Qi, ■ ■ ■ ,Q n +\) = (l,...,ra+ 1), 
(Si, . . . , 5 2rl +i) = (1, . . . , 2n+ 1) and multiply the result by a combinatorial factor (n + l)!(2n + 1)!. 
We first carry out the summations over gi , . . . , g 2 n an d fci , • • • , k n by following the analogous 
calculation for the ordered phase, see section 4.4. Finally, we carry out the sum over <72n+i using 
(D.4). This results in 



D 



(2n+2|2n+l|2n) 



h.C. = - 



2A 

! 



n 



Pi,...,p„GNS r=l 



A£ 4 

l-j- I 9 H (£-2te' Pr )[2ts' Pr 



x ^ J] sin(fc„ +1 ^) cos(2te fc „ +1 ) + . . . 



0<fc„ + iGNS 



4/ 4 

! - -- ^e H (e-2te' Pr )[2te' Pr -£} 



L L 



1 

Pi,...,p„SNS r=l 

x [£>(2|i|o) + -D(o|i|2)] +••• • (245) 
The sum over n can be again taken by inverting the steps used to express (120) in terms of (125), 
which gives 

oo 

^ £>(2n+2|2n+l|2n) + h.C. 
n=0 



[£(2|i|o)+£(o|i|2)]cxpf Yl ^U + K 2 ( P )(l-^(£ + e H (£-2te / p )[2te / p ~£})) 



- 0<pGNS 

[£> (2 |1| 0) +£>(0|1|2)] ns< s I s )ns ex P 



^ K 2 {p)(£ + e H (£-2te' p )[2te l p -£]) 



0<pGNS 



+ ... 



(246) 



In the last step we have only retained the 0(K 2 ) term in the exponent in order not to exceed the 
accuracy of our calculation. 
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4-5.6. Full Answer For the 2-Point Function in the Disordered Phase 

The above results suggests that the two-point function consists of two parts 

lim p xx {t,£) = F 1 (£) + F 2 (t,£). (247) 

L— 7-00 

The two kinds of contributions are 

• An oscillating time dependent contribution arising from the terms D( 2n +2\2n+i\2n) + h.c. 

F 2 {t,£)~2JA [ — sin(2te fe - ki) 

J-ir 1" £fc 

_ 2 f ^K 2 (p) (£ + 9 H (£- 2te' p )[2te' p -£})]. (248) 
Jo n 



x exp 



An exponentially small, time-independent contribution arising from the "diagonal" terms 

1 °° 1 ' 

^2 7T\ ^2 



NS (B\B) m ~nl 0<fc fc g 



x Ns(-k n ,k n ,...,-k 1 ,k 1 \a^ l+e a^ l \k 1 ,-k- L ...,k n ,-k n ) NS . (249) 

The 0(K ) and 0(K 2 ) have been evaluated in sections 4.5.1 and 4.5.3 respectively and are 
given by 

F 1 (£) « 2JA P p-— [l + 2K\q) At f ^K 2 (k) 
J-n 2tt e q [ J 2tt 

It is shown in Appendix E that the "pair ensemble" average (249) is equal to the average in 
the generalized Gibbs ensemble, which was previously calculated in [71]. Hence we conclude 
that 

9- 

1 + hh + v/(/i 2 - l){hl - 1) 



4(h 2 



"/" 

Jo 



F 1 (£)~C PP (£)exp 

= In (min[/io, hi]) — In 



2hh 



(251) 



where /ii = 
paper II. 



— 1+hha+ ^h+h a 1 ^ h ° ~ an< ^ ^ e ^ ar S c distance behaviour of Cpp(^) is determined in 



5. Scaling Limit of the Ising Model 

So far we have focussed on the quench dynamics in the transverse field Ising lattice model. In the 
vicinity of the quantum critical point at h = 1 a quantum field theory description applies, see e.g. 
[99]. As quantum quenches in integrable fields theories are of great current interest, we now present 
explicit expressions for the quench dynamics in the field theory limit. 

The scaling limit of the transverse field Ising chain is (a is the lattice spacing) [99] 

J^oo, h^r 1 , a ^ 0, (252) 
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while keeping fixed both the gap A and the velocity v 

2J|l-/i|=A, 2Ja = v. (253) 

In this limit the dispersion and Bogoliubov angle become 

e(q) = VA 2 + v 2 q 2 , (254) 
e h {q) arctan(^). (255) 

In our quench problem both the initial and the final magnetic field are scaled to the critical point, 
i.e. we need to take 

h Q -> 1 , J(l - h ) = A = fixed. (256) 
In this limit if-matrix turns into 

"arctan(^) - arctan(^) 



K{q) = tan 



Here the physical momentum is defined as 
k 



a 



-CO < q < CO. 



The Hamiltonian describing the scaling limit is expressed in terms of Majorana fermions as 

r °° dx 



H 



2n 



IV 



(tpd x ip - ipd x ip) - iAipip 



The order parameter in the scaling limit must be defined as 
a(x) oc (1 - h 2 Y^G^. 

where x — na a . It is customary to choose the normalization of the field a(x) such that 
km (0|a(sM0)|0) = -in 

X-S-0 | x |4 

which implies that 

a x ^ 2 1 /24 e l/8_ 4 -3/2 a l/8 (7(x) ( 

where 

A= 1.28242712910062... 



(257) 

(258) 

(259) 
(260) 
(261) 

(262) 
(263) 



5.1. Ordered Phase 

The result for the 2-point function after a quench within the ordered phase in the scaling limit is 



(Mt)HxM0)\Mt)) Kexp ' ' x 



(V'olV'o) 



f°°dq 
Jo n 



In 



l-K\q) 



1 + K 2 (q) 



[x9(2e , (q)t -x) + 2e'{q)t6{x - 2e'{q)t)]j , 

(264) 

where x — na n . This expression is well-defined because at large q we have 

A- A 



K(q) 



2vq 



(265) 
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In the stationary state we have 

a(:> 

t^oo (V'olV'o) ~ ^ V (J ' 

where 



1 = A + Ao_^_ 

s 2u v v ' 



5.2. Disordered Phase 

For quenches within the disordered phase we can take the scaling limits of the results (248) and 
(249), which give 

(Mt)\<r(x)a(0)\Mt)) 



(V'olV'o) a + T2 ^ ,X ^ 1 

where the time-dependent part dominates except at very late times and is given by 



(268) 



F 2 (t,x) « 2v f 

J — ( 



x cxp 



dq K(q) 
2tt e(q) 
dq 



sin (2te(q) — qx) 
l-K\q) 



In 



l + K 2 {q) 



The stationary state component is 



T\ (x) oc exp 



[x6{2s'{q)t - x) + 2e'{q)te{x - 2e'{q)t)] ) . (269) 

(270) 



where now 



_ A + A _ VAA + min(A , VAA ) 

2v v v 



(271) 



6. Conclusions 



In this work we have derived analytic expressions for the time evolution of one and two point 
functions in the transverse field Ising chain after a sudden quench of the magnetic field. To do so 
we have developed two novel methods based on determinants and form factor sums respectively. The 
former is applicable to quenches in models with free fermionic spectrum and our analysis generalizes 
straightforwardly e.g. to the spin-1/2 XY chain in a magnetic field [88]. Results obtained by this 
method are exact for asymptotically large times and distances in what we call the space-time-scaling 
limit (t, I — > oo keeping their ratio fixed). The form factor approach is applicable more generally to 
integrable quenches [32] in integrable quantum field theories [100] such as the sine-Gordon model. 
It is furthermore straightforwardly extended to the study of non-equal time correlation functions 
[101]. The form factor method provides approximate results that become exact in the limit of 
small quenches, defined by the requirement that the density of excitations (of the post-quench 
Hamiltonian) in the initial state is low. We observe that the difference of the form factor and exact 
results for quenches within cither the ferromagnetic or paramagnetic phase are generally very small, 
except for quenches originating or terminating in the close vicinity of the quantum critical point. 
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Appendix A. Diagonalization of the Transverse Field Ising Model 

In this Appendix we summarize the diagonalization of the TFIM [102] with periodic boundary 
conditions 

L 

H(h) = -jjy*°hi + h < j ji (at) 

3 = 1 

where I is even, cr" are the Pauli matrices at site j and 

=of , a = x,y,z. (A.2) 

The dimensionless constant h describes the coupling with an external magnetic field Jh. The 
quantum Ising chain is mapped to a model of spinlcss fermions by means of a Jordan- Wigncr 
transformation. Defining — [af ±z<rjf)/2 we construct spinless fermion creation and annihilation 
operators by 



4 = 
-i 

3 = 1 



The inverse transformation is 

j-i 

o*i = 1 - 2c} Cj , <r? - 1] (1 - 2c] Cl ) ( Cj + cj) . (A.4) 
l=i 

The Hamiltonian can be expressed in terms of the fermions as 

L-l L 

H {h) = - j J2^4 - c 3ii c 3+i + c ]+i\ - Jh J2 c ^4 ~ 4 C 3 

3=1 3=1 

- Je™*(c L - c{)( Cl +c\), (A.5) 

where 

L 

N = J2 c h- ( A - 6 ) 

3 = 1 

As [H, e 47rAr ] = we may diagonalize the two operators simultaneously. The Hamiltonian is block 
diagonal H = H e H D , where H e / a act on the subspaces of the Fock space with an even/odd 
number of fermions respectively. 
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Appendix A.l. Even Fermion Number 

In the sector with an even number of fermions we have e l7rN = 1 and the Hamiltonian can be 
written in the form 

L L 
H e (h) = - J^[ct - cj][c j+1 + ct +1 ] - JhY,c 3 c] - c]c 3 , (A.7) 

where we have imposed antiperiodic boundary conditions on the fermions 

cl+1 = -ci. (A.8) 
The Hamiltonian H e is diagonalized by going to Fourier space 

1 L 

c(k n )= £ Cj . e **»i, ( A - 9 ) 

where fc„ are quantized according to (A.8) 

,„ = M=±M, „ = -£,...£-!. (A,0) 

The antiperiodic sector is commonly referred to as Neveu-Schwarz (NS) sector. Following this 
nomenclature we introduce the notation k £ NS to describe the set (A. 10). Introducing Bogoliubov 
fermions by 

c(k n ) = cos(0 fen /2)a: fen +ism(6 kn /2)a*_ kn , 

cH-k n ) = ism(6 kn /2)a kn + cos(0 kn /2)al kn , (A.ll) 
where the Bogoliubov angle fulfils 

h — p ik 

e i9 » = , - (A.i2) 

^1 + h 2 -2hcosk 
the Hamiltonian becomes diagonal 



-l 



H e (h)= £ ( fc «) 



(A.13) 



Here the dispersion relation is 

e fe = 2Jy/l + h 2 -2hcos{k). (A.14) 
A basis for the Fock space in the sector with even fermion number is then given by 

2m 

|fci,...,fc 2m ;/i)NS = II a U 0; ^ NS ' fc i eNS ' ( A - 15 ) 

where the fermion vacuum |0; /i)ns is the state annihilated by all a kj (j = — f, — 1). 
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In the sector with an odd number of fermions we have e l7rN = — 1. The Hamiltonian can again be 
written in the form 

L L 

H a (h) = - Jj>j - Cj }[c j+1 + ct +1 ] - Jh^jc] - c]c 3 , (A.16) 
but now we have to impose periodic boundary conditions on the fermions 



cl+1 = ci- 
la Fourier space we therefore now have 

L 



Cj e 



(A.17) 
(A.18) 



2irn 

Pn = , n - 



1. 



(A.19) 



where p n are quantized according to (A.17) 

L L 

L ' " 2' "2 

The periodic sector is known as Ramond sector and we will denote the set (A.19) by p n G R. 
Defining Bogoliubov fermions a Pn for p n ^ by 

c(p n ) = cos(9 Pn /2)a Pn +ism(9 Pn /2)a<_ Pn , 

J{-p n ) = isin{6 Pn /2)a Pn + cos(6» Pri /2)al pn , 

we can express the Hamiltonian as 



(A.20) 



H (h)= ]T e{ Pn ) 



t 1 



2J(1 - h) 



"=~2 



A basis of the subspace of the Fock space with odd fermion numbers is then given by 

2m+l 

\pi, ■ ■ ■ ,P2m+i;h) = Yl a\ 3 10; h) R , p 3 E R, 



(A.21) 



(A.22) 



3=1 

where the fermion vacuum |0)r is the state annihilated by all a Pj (j = — ^, . . . , \ — 1). 

Appendix A. 3. Paramagnetic Phase h > 1 
Here the ground state is 
|0)ns- 

A complete set of states is then given by 

2m+l 

\pi,---,P2m+i;h) R = Yl al.\0; h) R , 

2m 



(A.23) 



|fci, . . . , k 2m ; h) NS = Yl a lj l 0; ^ NS • 

Pj GNS 



(A.24) 
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The Hamiltonians can be written as 

H e (h) = £ ( fc «) a L a ^ +Eo S (h) , 

fc„eNS 

H (h) = J2 <Pn) + E oW> (A.25) 

where E-(h) = -\ E g£a e(<z), a = R, NS. 
Appendix A. 4- Ferromagnetic Phase h < 1 

As the zero momentum mode has negative energy it is useful to perform a particle-hole 
transformation 

a — ► aj . (A.26) 

Redefining the Ramond vacuum as the state that is annihilated by all a Pn after the particle-hole 
transformation we can construct a complete set of states as 

2m 

\k 1 ,...,k 2m ;h) R = U<4jO;/i) R , 

2m 

|pi,...,P2m;fe>NS= J] "pJ ; fe >NS ■ (A.27) 

The Hamiltonians are then again given by (A.25). For large L we have EQ S (h) — E^(h) = C(L _1 ), 
so that there are two low-energy states 

\0;h) K , |0;/i)ns. (A.28) 

As long as L is finite the ground state is |0;/i)ns- O n the other hand, in the thermodynamic 
limit the states (A.28) become degenerate and by spontaneous symmetry breaking one of the two 
combinations 

-^[|0;fc) R ±|0;/i>Ns] (A.29) 
is selected as the ground state. 



Appendix B. Initial State 

As described in Appendix A the Hamiltonian H(ho) can be diagonalized by a Bogoliubov 
transformation. Let us denote the corresponding Bogoliubov fermions by &k, the Bogoliubov angle 
by 9® and the NS vacuum by |0;/io)ns- Similarly the Hamiltonian H(h) is diagonalized by the 
Bogoliubov transformation (A. 11) and its lowest energy state in the even fermion sector is |0; /i)ns- 
As both sets of Bogoliubov fermions are given in terms of the same spinless fermions Cj and ct 
(j = 1, . . . , L), they can be expressed in terms of one another by 

a fc„=cos^ — -)a kn +i sin — -ja_ kn . (B.l) 
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As both sets of fermions can be used to construct a basis of the even Fock space we can express 
|0;/iq)ns in the form 



INS 



n=0/c(i),...,fc(n)£NS j=l 

Using the expression (B.l) in the condition 

"kjO; h ) m = 0, 

allows determination of the coefficients f k m ... fc(«)- A simple calculation gives 



\0;h ) 



NS 



Aft 



NS 



exp 



PENS 



|0; /i)ns, 



where A/ns is a normalization constant and the function K (k) is given by 

'0 k -e^ 



K{k) = tan( 

The equations of motion for a k imply that 

a k (t) = e U£k a k (0), 

so that [61] 

c -itH e (h)i n . h \ |-B(*))nS 

VNS(i*|-D) NS 



where 



\B(t)) m = e~ tm o exp 



0<pGNS 



|0;/i) NS - 



Similarly one can show that 



e-* g °WlO;fto) R = j£Mg , 
^{B\B) R 



where 



\B(t)) R = e- UE o exp 



0<pSR 



\0;h) R . 



(B.2) 
(B.3) 
(B.4) 

(B.5) 
(B.6) 

(B.7) 
(B.8) 
(B.9) 
(B.10) 



The states (B.8) and (B.10) are of the same form as boundary states in integrable scattering theories 
[94, 9]. 

A physical interpretation of the function K (k) (B.5) is obtained as follows. The density of 
post-quench Bogoliubov fermions a\a k in the initial state is given by 



(0; h \a\a k \Q- h ) a = a (0; h\ ]J (1 - iK{k)a k a- k ) a\a k (l + iK(k)al k a\}j [0; h) 



fcGa 



K 2 {k) 
1 + K 2 (k) 



, a = NS,R. 



(B.ll) 
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This implies that in the case where K(k) is uniformly small in k we have 

(*ol4"fel*o> - K 2 (k) + O (K 4 ) , (B.12) 

where |\&o) is the initial state of our quantum quench. Physically the small parameter characterizing 
the expansion in powers of K(k) is therefore the density of excitations of the post-quench 
Hamiltonian H{h) induced by the quantum quench. 

Appendix C. Proof of the product formula 

We need to evaluate the product of the 2x2 matrices 

2n 

n 2n ({a t }) = Y[M{ ai ), with M(a) - Aa x + Ba v e iaa * . (C.l) 
»=i 

Compared to Eq. (85) in the main text we have A — n x (k ), B = \n±(ko)\, <Xj = 2ei_it and 
we choose (without loss of generality) n± = y. Using the algebra of the Pauli matrices, it is 
straightforward to calculate 

n 2 (oi, a 2 ) = M(ai) • M(a 2 ) = A 2 l + B 2 e l< - a2 - ai ^ + iABa z {e ta ^ - e iaia *) , (C.2) 

where I is the 2 by 2 identity matrix. A slightly longer exercize is required to calculate P±4 and to 
obtain 

4 

n 4 = ^ 4 I + L4 3 Bo- z ^(-l)V a * CTa = + A 2 B 2 {-l) jl+j2+1 e i ^- a »^ (C.3) 

j=l l<h<32<i 
+ iAB 3 a z ^ — iyi+h+j3 + l e i(a j3 -a j2 +a jl )rT x _|_ Qi e i{a i -a 3 +a 2 -a 1 )(T x _ 

1<J1<J2<J3<4 

From these two first examples, it should be clear that the general structure of n 2n is 

2n p 

n 2n = Y, A2n ~ PBP ^ P E (-l) ELlJfc exp(^(-l) p -^*)- (C4) 

p=0 l<j 1 <j 2 <...j p <2n fc=l 

Having this conjecture, it is straightforward (but require some algebra) to prove it by induction 
showing that it is compatible with the recurrence relation 

ri2n+2(ai, • • • «2n+2) = n 2 „(ai, . . . a 2n ) ■ n 2 (a 2 „+i, a 2n+2 ) • (C5) 
Appendix D. Useful relations 

Lemma 1: For any function f(z) that is 2tt periodic and analytic in a strip around the real axis we 
have for < k € R 



/(gn)e 2« te(g „) 

L qn ^ S >0 (cOSfc-CO S(?n )2 



\-te\k) 



e 2lt£k f(k) + ^e 2te * t f'(k) 



sin 2 k 2sin 2 (fc) 
dz f{z)e 2l£ ^ t 



I± (Dl) 

J 27r(cosz-cosfc) 2 (l + e liz )' v 

where the integration is along a closed contour encircling the interval [0, 7r]. 
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Lemma 2a: For any function f(z) that is 2tt periodic and analytic in a strip around the real axis 
we have for < fc £ R 



2ite(q„) 



1 



L (cosfc — cos g„) 2 



6i cot k 



te'(fc) 



f(k) 



sin 2 k 



„2ie k t el 



+ i- 



/'(*) 



2sin (fc) 



/(*)e 2fe M* 



27r (cos z — cos fc) 2 (l — e lLz ) 



0(L- 



(D.2) 



The analogous equation for momenta in the NS sector is 
Lemma 2b: For any function f(z) that is 2ir periodic and analytic in a strip around the real axis 
we have for < k £ NS 



1 

L 



E 



f{qn)e 



2ite(q„) 



,„eNS 



>o (cosfc - cos q n ) 2 



L 
12 



1 - 



6i cot k 



te'(k) 



f(k) 



sin 2 fc 



f(z)e* 



2ie(z)t 



0(L _1 ). (D.3) 



2sin 2 (fc) / 2tt (cosz - cosfc) 2 (l + e lLz ) 

Lemma 3: For large j 3> 1 and any function f(z) = f(z + 2w) that is analytic in a strip around the 
real axis we have (fc £ R) 



Um 1 E /fa,)e«- _ 

L^co L COS fc — COS Q n 

where 7 is a positive constant. 

Proof: Using contour integration we have 

1 \ - f(q n )e^" f dz 



-ikj 



f(-k) 



a ikj 



f(k) 



2i sin fc 



+ 



(D.4) 



L ^ cos fc — cos q n 



2n (cos z — cos fc)(l + e lLz ) 



-ikj 



+ 



f(-k) 



f(k) 



2 ,sinfc (°- 5 ) 
where we have used the quantization conditions (A. 10), (A. 19) and where the integration is along a 
counterclockwise countour encircling the interval [— 7r, tt] in a manner such that no singularities of 
f(z) lie within it. The contribution of the part of the contour below the real axis tends to zero in 
the limit L — > 00 because |e lz ( L±:, '| >• 1. As j > the part of the contour above the real axis can 
be deformed as shown in Fig. Dl. The contributions of the pieces parallel to the imaginary axis 



-JC 



7T 



Figure Dl. Deformed integration contour. 
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cancel due to the periodicity of f(z), while the remaining part has the exponentially small bound 
given in (D.4). 

Lemma 4-' For any function f(z) that is 2ir periodic and analytic in a strip around the real axis we 
have 



E 



/fan) 



1 



L qf^iS ( COsk - COS 1n) 2 



e-ikifi-k) + e ik 'f'(k) 



+ 



4 2 

Proof: Using contour integration we have 



2i sin 2 k 



f(k) 



f(q n )e 



1 y 

L qt&is (cosk-cosq n ) 2 



-ikj 



f(-k) 
sin 2 (A;) 



e ikj f(k) 



sin 2 k 



4 2 



+ o(e~ 



13 



(D.6) 



+ i 



e- ik *f'(-k) + e ik 'f'{k) f dz 



2sin 2 (fc) 



f{z)e^ 



2n (cosz — cosfc) 2 (l + e tLz ) ' 



(D.7) 



where the integration is along a closed contour encircling the interval [— tt,it]. As for lm(z) < we 
have 



lim e i{j - L ^ 

L— >oo 



0, 



(D.. 



only the upper part of the contour contributes in the L — > oo limit. Using that f(z) is analytic 
in some strip around the real axis the part of the contour above the real axis can be deformed as 
shown in Fig. Dl. The contributions of the pieces parallel to the imaginary axis cancel due to the 
periodicity of f(z), while the remaining part is exponentially small in j. In the limit of large L and 
j the first term in (D.6) can be neglected, so that 



E 



f(q n )e 



-ikj 



f(-k) 



Jkj 



f(k) 



(cos A: — cosg„) 2 



sin (A) 



4 2 



(D.9) 



Lemma 5: For any function f(z) that is 2tt periodic and analytic in a strip around the real axis we 
have 



1 



E 



hm 

L->oo L £ -~' cos k — cos q 



f(q) e 24<Te («)* iaf(k) e 2l<T£fct r +al ° dz f(z) e 2lae ^ t 



2 sin k 



aid 



2tt cos k — cos z 



(D.10) 



where a = ± and < k e R. 



Appendix E. "Pair Ensemble" Averages 



The "diagonal" terms in the Lehmann representation based on the eigenstates of the post 
quench Hamiltonian H(h) give rise to a time-independent contribution to the expectation value 
( x f , o(0)|0| v I'o(0)), which we call "pair ensemble". It is defined as the average 



1 ^ 2 i 

^ PE = Ns(B\B) m E n\ E 

ins\ i ' n= o o<fei,...,fc„eNS 



x ns ( k n ,k r , 



U k2 ^ 

3 = 1 

. , -hi, fci; h\(D\ki 7 -fci . . . , k n ,-k n ; h) NS . 



(E.l) 
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Averages of the form (E.l) can be represented using a density matrix as (O)pe = tr(ppE0), where 
p PE = (l-n )(l-n 7r ) j[ ( + n k nA. (E.2) 

CKfceNS 1 + ^ W 

Here n q = a q a q are the Bogoliubov fermion number operators in momentum space and I k is the 
identity in the subspace with momentum k. The only non- vanishing expectation values of the pair 
ensemble are 

(n qi n g2 . . . V )pe = II ifSh ' (E - 3) 

\ qj \es 

where S is the set of all momenta qj with mutually distinct magnitudes, i.e. 

< |<?r| + \q*\ < n Vq r ,q s e S. (E.4) 

We note that (n k n- k O) = (n k O) since the ensemble describes pairs of particles with opposite 
momenta. 

As shown in section 3 lattice spin operators can be expressed as products of the Majorana 
fermions aj, a v -, which are related to the Bogoliubov fermions diagonalizing the post-quench 
Hamiltonian H{h) by 



ikj p i6 k /2 

3 Vl ' 



k 



aV o = 4= E e- ikj e~ i9k/2 i \al k - a k ] . (E.5) 

As the only non-zero expectation values in the pair ensemble are of the form (E.3) we can 
express a general average in the form 

2n ^ n s 

IK;) PE = ^ E £4 s) ( fc i'---'^)<n^>E' (e.6) 

r— 1 ki,....k n s— it— 1 

•e X^(^i) ■ • • > &n|i) are well-behaved functions of the momenta. 
The generalized Gibbs ensemble for the Ising model is defined as 

(0) GGE = — l—trfe" ©I , (E.7) 

^GGE L J 



where Zqge = tr 



and 

A, = - log (K 2 (q)). (E.8) 
By taking the trace over a basis of eigenstates of H(h) we conclude that the only non-zero averages 



are 

-2 



\ n qi n q2 ■ ■ ■ n q m ) GGE = [[ 1+R 2^ ' 



(E.9) 



The corresponding density matrix is 



_ „ l k -(l-K^(k))n k 
PGGE 1 + ^ 2(fc) • 
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We observe that (n^nj-j • • • rik m ) GGE = (it^QGE ( n k2 ' ' ' n fe m )cGE f° r an y se ^ 01 distinct momenta 
fci, ^2, • • ■ , k m , and unlike in the pair ensemble Wick's theorem applies in the generalized Gibbs 
ensemble. General averages can be expressed as 

2n 1 n s 

(li<X G E=L^ E £4 S W--,Mi)<n^>GGE, (E.ll) 
r—1 k±,...,k n s—Q u—1 

where the functions X^\k\, . . . , k n \j) are the same as in (E.6). As the averages (n«=i ^Ogge 
and (n«=i n fc u )pE are the same unless at least two of the \k u \'s coincide, and such contributions 
are suppressed by factors of 1/L, we conclude that 

2n 2n 

lim (T]af) = lim (T\ ab f) ■ (E.12) 

r—1 r—1 

The above arguments show that averages of operators that are local in space are the same in both 
ensembles. However, non-local operators have in general different averages, e.g. (nfcn_fc) GGE = 
(n k n- k )-p E . 
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